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Methods in Transient Stress Analysis 


R. L. BISPLINGHOFF,* G. ISAKSON,?7 anp T. H. H. PIANt 
Massachusetts Institute of Technology 


SUMMARY 


A general method of approach to the problem of determining 
transient stresses in damped and undamped elastic structures is 
outlined, and methods for the approximate treatment of contin- 
uous elastic structures are discussed. The formulation of the 
problem is based on Lagrange’s equation of motion and the use 
of displacements of natural modes of vibration as generalized 
Two methods are found suitable for solving the 
The first is an 


coordinates. 
equations of motion of the transient problem. 
application of the Laplace transformation which has been ex- 
tensively used in the solution of transients in electrical systems. 
The second employs the principle that the transient response 
to an arbitrary forcing function can be expressed in terms of the 
fesponse to sinusoidal forcing functions. The response due to a 
sinusoidal forcing function is defined as the mechanical admit- 
tance. Two methods of approach for determining the stresses 
corresponding to the deformation in an elastic structure are 
outlined. In the first, the stress is considered to be a super- 
position of the stresses corresponding to displacements of each of 
the normal modes. The stresses in each mode are calculated by 
Using the inertia loadings associated with the natural vibrations. 
In the second, the total transient stress is considered as the sum 
of the stresses that would result if the system were restrained 
against vibrating and the additional stresses that result from the 
vibratory motion. 


(1) INTRODUCTION 


I" THE DETERMINATION of the design loads on an air- 
plane structure in an accelerated condition, it is 
usually assumed that the airplane is perfectly rigid. 
Structural components designed by loads computed on 
this basis have been known to fail in service. It 
has been established in some of these cases that 
failure was due to a dynamic overstress ignored by the 
designer. 

This dynamic overstress is a consequence of the lack 
of rigidity of the structure. External loads that are 
fapidly applied not only cause translation and rotation 


of the airplane as a whole but also tend to excite vibra- 
a 
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tions of the structure. It is the additional inertia 
forces associated with these vibrations that account for 
the dynamic overstress. 

In analyzing the airplane as an elastically deformable 
body, the pseudo-static treatment previously employed 
is no longer sufficient. The problem becomes one of 
studying the transient motion and deformation of the 
airplane and determiniag the resulting variation of 
stress with time. 

A number of papers have been written on this sub- 
ject. It is usually assumed, however, that the reader is 
familiar with the mathematical concepts and methods 
that are involved. This information must, in fact, be 
culled from a variety of sources, where it exists in a 
connection generally somewhat removed from that of 
the present application. It therefore appeared de- 
sirable to the authors to make a clear statement of the 
problem and to present a systematic introduction of the 
mathematical concepts and methods that are required 
foritssolution. This is attempted in the present paper, 
and some additional methods, new to the present ap- 


plication, are introduced. 


DESCRIPTION OF THE SPACE CONFIGURATION OF 
AN UNRESTRAINED ELASTIC SYSTEM 


(II) 


In the treatment of problems in dynamics, the proper 
choice of coordinates is important, since the labor re- 
quired to obtain a solution depends largely upon this 
If the coordinates are chosen so that they are 
completely independent—that is, they are not related 
by equations of constraint—the analysis will generally 
be greatly simplified. Such coordinates are known as 
generalized coordinates. In a system described in this 
way, the number of degrees of freedom is equal to the 
number of generalized coo dinates required. 

A rigid body moving freely in space possesses six de- 
grees of freedom, and its position may be established by 
three linear displacements and three angular displace- 
Its motion may be determined by the solution 


choice. 


ments. 
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of total differential equations in which time is the only 
independent variable. 

A continuous elastic body, on the other hand, 
possesses an infinity of degrees of freedom, since it is 
not possible to represent a general deformation of the 
body in terms of a finite number of coordinates. The 
motion of such a body or system may be treated an- 
alytically by means of partial differential equations or 
integral equations, which contain unknown functions 
of both time and space coordinates. Except in the 
simplest cases, however, these equations are extremely 
difficult to solve exactly. In order to make such a 
problem tractable, it is necessary to find a suitable 
system of generalized coordinates which will make it 
possible to formulate the problem in terms of total 
differential equations. Two different approaches are 
often used in accomplishing this. 

In the first, the elastic body is regarded as a system of 
discrete mass particles that are elastically connected. 
In an exact treatment there would be an infinity of 
such particles. It is often satisfactory, however, to 
approximate the actual system by lumping portions of 
the mass of the body into rigid bodies that are then 
considered to be connected by weightless elastic con- 
nectors. The elastic properties of this system may be 
expressed in the form of a set of influence coefficients 
linking every pair of masses. If the system is free in 
space, it is usually convenient to locate the system as a 
whole by specifying the position in space of a set of 
orthogonal axes which is fixed in some way to the system 
(for example, it may be fixed to one of the masses) and 
locating the individual parts with respect to these axes. 

In the second approach, any deformation shape of 
the continuous elastic body is considered as a super- 
position of certain assumed deformation shapes that 
are compatible with the boundary conditions. Each 
assumed deformation shape or mode represents a de- 
gree of freedom, and the multiplier that determines 
the amount of its contribution to any general deforma- 
tion represents the generalized coordinate correspond- 
ing to that mode. In an exact treatment an infinity of 
such modes would be required. In many practical 
problems, however, a relatively small number of modes 
will satisfactorily approximate the significant deforma- 
tions of the system if the modes are carefully chosen. 
In addition to the deformations, there may, of course, 
be displacements of the system as a whole, which may 
again be represented by locating axes that are fixed in 
some manner to the system. 


(III) THE Equations oF MoTION 


By describing the configuration of the elastic system 
in terms of a finite set of generalized coordinates selected 
in either of the two ways described above, it is possible 
to write the equations of motion of the system in the 
form of total differential equations, with time as the 
sole independent variable. This can be done by ap- 


plication of D’Alembert’s principle and the principle of 
virtual work through the form known as Lagrange’s 
equations." These equations are written as folloys 
with each equation corresponding to a generalize 
coordinate of the system: 
d (=) oT , OV 
dt \Oh, Oh, dh, 
where h, is the displacement in the rth generalized co. 
ordinate, h, is the first derivative of h, with respect to 
time, 7 is the kinetic energy of the system, JV’ is the 
potential energy of the system, and //, is the general. 
ized force corresponding to the generalized coordinate }, 
and is the work done by all forces other than inertia or 
conservative forces per unit of virtual displacement jy 
the rth coordinate. 

IT, can be defined mathematically as follows: 

— él, . 
bh, 

where 6/V, is the work done by all forces other than 
inertia or conservative forces in a virtual displacement 
6h,. 

The kinetic energy of the system may be determined 
from the fundamental formula, 


T = (1/2) / v2 dm (3 
M 


where v is the linear velocity of the element of mass 
dm and the integration is carried out over the whole 
mass of the system. 

The following assumptions, which will ensure that 
the equations are linear, are valid for the particular 
problems with which we are concerned: (a) The de- 
formations of the system are small, and (b) angular 
displacements of the system as a whole are small. 

With the introduction of these assumptions, the 
kinetic energy becomes a homogeneous quadrati 
function of the generalized velocities h, . . . h,, of the 
form,” 


T= ee mM jh ih; 4 
2 i=l j=! 


The quantities m,; are inertia coefficients and are con 
stants of the system. They exhibit the property mi; = 
M ji. 

The potential energy V is a measure of the ability 
of the conservative forces to do work upon the system. 
In the problems with which we are concerned, thest 
forces are elastic forces, and the potential energy is the 
strain energy associated with the deformation of the 
structure. Gravity forces do not, in general, enter the 
problem, since they are steady forces and, consequently, 
make no contribution to the disturbed motion. As 4 
consequence of the assumption of small deformations, 
the potential energy can be expressed as a homogeneous 
quadratic function of the generalized displacements.” 
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METHODS IN 


n n 


V = (1/2) } >> Rahs; (5) 
#=1 j=l 

where k,; are elastic constants of the system and ex- 

hibit the property ki; = kj, as a consequence of Max- 

k,; will be 

zero When 7 or j refers to a coordinate for the system as 


well’s theorem of reciprocal deflections. 


a rigid body, since displacement in such a coordinate 
does not involve deformation of the system. 

Introducing expressions (4) and (5) into Eq. (1), the 
equations of motion may be written, 


> mryhjy + > kak, = H, (¢¥ = 1,2,3...) (6) 
j=l j=1 


where i; is the second derivative of 4; with respect to 
time. 


THE EQUATIONS OF MOTION IN TERMS OF 
NORMAL COORDINATES 


(IV) 


The equations, as they stand, contain inertial and 


elastic coupling terms. It would however, 
that their solution would be simplified if these coupling 


It is, in fact, possible to de- 


appear, 


terms were not present. 
termine a system of coordinates such that there is no 
inertial or elastic coupling. Such coordinates are 
known as normal coordinates. 

If the expressions for kinetic and potential energy are 
of the form given in Eqs. (4) and (5), it is always pos- 
sible to find a linear transformation of the coordinates 
such that the kinetic and potential energies, when ex- 
pressed in terms of the new coordinates q,, reduce to the 
form? 

n 


T = (1/2) >> M.q,? (7) 


, l 


n 
VY = (1/2) z. M,w,*q,? (S) 

7 1 
where the quantities .1/, and w, are constants of the 
system and are defined below. This transformation can 


be represented by the equations 
n 
h, = > ae, ( -@ 1,2; 35..8) (9) 


It can be shown? that the quantities w, and \;”’ may 
be determined by substituting the solution 


h, = \,°" sin (w,t + V,) 


into the differential Eqs. (6), made homogeneous by 
setting the right-hand side equal to zero and solving the 
associated characteristic equation. This yields a set of 
n characteristic roots w,, and associated with each there 
isa set of n characteristic values \,"” (i = 1, 2,... 7). 
Physically, w, is the circular frequency of the rth prin- 
cipal oscillation and is a natural frequency of vibration 
The n quantities \;"" associated with 
#, are the mode coefficients determining the configura- 


of the system. 
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tion of the system in the rth principal oscillation or the 
rth natural mode of vibration. 

There will, in the case of an unrestrained body, be 
normal coordinates locating the principal inertia axes 
These represent modes of zero 
frequency involving no deformation of the body. The 
vibration modes, on the other hand, involve no dis- 


of the body in space. 


placement of the principal inertia axes but only defor- 
mation with respect to these axes. 

When the body is replaced by a system of lumped 
masses, the initially chosen coordinates /; represent dis- 
placements of these masses, and the mode coefficients 
\,’ determine the relative displacement of the masses 
corresponding to each normal mode. On the other 
hand, when the body is considered to be continuous and 
h, represent displacements of assumed deformation 
determine the relative con- 


(Y) 


modes, the coefficients X, 
tributions of these modes to the normal modes, which 
may then be represented as continuous functions of the 
space coordinates designated by ¢"”. 

The determination of the characteristics of the prin- 
cipal oscillations of a continuous system can also be 
considered from the point of view of the solution of a 
homogeneous linear integral equation. In Appendix 
(A) an integral equation is derived for the case of an 
unrestrained elastic airplane, and it is shown that there 
are associated with this equation an infinite number of 
characteristic values w, and characteristic functions ¢”. 

The quantity J/,, called the generalized mass cor- 
responding to the rth normal coordinate, is now defined 


as follows: 


M,= > Da mais? (10) 
s=1 9 


where mm; are the inertia constants of the system intro- 
duced earlier. 
It can be shown® that the m sets of mode coefficients 
satisfy the orthogonality condition 
n n 
y~ DY marMra® =0 (r # s) (11) 
i=1 j=1 
In the case of a system of lumped mass particles, 
miy = 0; mj, is the magnitude of the 7th mass, which 
may be written simply m,; and we have, 


nN 


M,= > ma,” (12) 
i=1 
while the orthogonality condition becomes 
n 
y mAOru” =0 (r #5) (13) 


t=1 


Similarly, in the case of a continuous system, the gen- 
eralized mass may be written 


M, = | od dm 
ul 


and the orthogonality condition is 


(14) 
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J obdm = (0) (r y= §) (15) 
M 


The integration in each case is carried out over the 
whole mass of the system. 

The mode coefficients \,;“? corresponding to each 
mode are not determined absolutely in the analysis; 
only their relative values may be determined. The 
process of determining absolute values for the coeffi- 
cients is arbitrary and may be carried out in a number 
of different ways. We may, for example, assign an 
arbitrary value, perhaps unity, to one of the A;°”’s 
and determine corresponding values for the others, or 
we may assign some arbitrary value to /,. This proc- 
ess, known as normalizing the modes, will determine 
the absolute values of qg,, the displacements of the nor- 
mal modes. It is, in effect, a matter of deciding what 
shall constitute a unit displacement of each mode, since 
the displacement of a normal mode represents a con- 
figuration of the whole system and not merely the 
simple displacement of a single point. 

It is often convenient to normalize the modes so that 
sach generalized mass MM, is equal to M, the total mass 


of the system. ‘ 
n n 
M= > Dd maa (16) 
é=1 j=l 


If it is assumed that the modes had previously been 
normalized in some other way, with coefficients A,“ 
describing the mode shapes, we have 


n 


n 
M,= > my; AyOA, 


s=1j=1 


(17) 


There is a factor, which we designate K,, that relates 
the \,"”’s to the A,’s, so that 


4” = KA (18) 


Substituting this into Eq. (16) and solving for K,, we 
obtain 


(19) 


n n 
! # = Vt ys >» mA {PAs 

i=1j=1 
In the continuous mass representation, this may be 
written 


K, = VM/fydm (20) 
where y is the function describing the normal mode 
shape originally, @ is the final normalized mode shape, 
and @ = K,y™. 

It is seen that K, depends only upon the distribution 
of mass and not upon the absolute mass of the system. 
The advantage of this type of normalization is that only 
one dimensional mass parameter appears in the equa- 
tions of motion, and it does so explicitly. All other 
quantities depend upon the distribution of mass and 
not upon its absolute value. 
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Substituting Eqs. (7) and (8) into Lagrange’s equa- 
tions [Eqs. (1)], we obtain the equations of motion ix 
terms of normal coordinates, 


Min = OQ; 
M ylin = On 
Pe 9 a 
M,, + 19m+1 i Mn 41%m + 1Qm4+1 = Om 1] az 
Mudn + Myeon?Qn = Q, 


Ym locate the system as a 
- Qn deter- 
mine the displacements of the deformation modes. 


where the coordinates q. . . 
whole in space and the coordinates g,, +; 


(V) SOLUTION OF THE EQUATIONS OF MOTION jw 


TERMS OF NORMAL COORDINATES 


The method to be employed in the solution of the 
equations of motion depends somewhat on the nature 
of the generalized force. When this force derives en- 
tirely from arbitrary external forces of known time 
history so that it is independent of the motion of the 
system, the equations are all independent of each other 
and each can be solved separately. 

In this simple case, with initial conditions q,(0) = 
g-(0) = 0, the solution may be obtained by first deter- 
mining the response to a unit-step function and then 
substituting the result into Duhamel’s integral, to 
yield? 


1 t 
2 = : ), sin w,(¢ — (22 
q(t) med | Q,(r) sin w,(t t)dr 


When the problem is such that the generalized forces 
derive not only from arbitrary external forces but also 
from forces that depend upon the motion of the sys- 
tem, the solution is more complex. The equations 
may then contain terms coupling the coordinates, so 
that a simultaneous solution of the set of equations may 
be necessary. 

Two general methods of approach for such problems 
are considered. The first involves the application of 
the Laplace transformation; the second makes use of 
the concept of mechanical admittances. 


(1) Solution by the Laplace Transformation 
The generalized force Q,(t) can usually be put in the 
form 


~ Gn) + 
O;?(t) (23 


Q(t) _ OM (1 7+ Qn op eee Gng M1 . 


where Q," represents that part of the generalized force 
which is due to the motion of the system and Q;"(/) Is 
the generalized force corresponding to the disturbing 
force F(t) which is responsible for the transient motion. 
It is assumed here that the applied forces do not intro- 


duce nonlinearity into the problem. 
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The equations of motion are 


ie) 
(24) 


. Ins Vi cee Gas qi . 


Mai — Mwirdi > Od ates 
G=_.25...0 


0,"(t) 
The Laplace transformation® of a function f(t) is an 
operation defined as follows: 


eff} = f(r) = f e-P" f(t) dt 
0 


where ef } denotes the operation of transformation. 
This yields a function, f(p), of the variable p, which is 
ysually called the Laplacian operator. 


(25) 


By applying this transformation to all the terms of a 
linear differential equation, we obtain an algebraic 
equation. Applying an inverse transformation to the 
solution of this algebraic equation, we obtain in cor- 
rect form the solution to the original differential equa- 
tion. The method is particularly useful when a num- 
ber of equations must be solved simultaneously, and 
more conventional methods of solution are not prac- 
tical. 

Applying the Laplace transformation to Eq. (24), we 
obtain 


Mpa(p) + Mwegi(p) = LL0(qM.-- de 
Gi. - + Qn Gr. ++ Gn) + OP(O} 


(26) 
where 
gi(b) = Liqilh)s 


As a consequence of the assumption of linearity of 
the applied forces, Eq. (26) can be put in the form 


(pe? + w2)Mi — Qi (p) bap) — 
D Qu" (P)aNb) = Q(P) (27) 
jF1t 
where Q;;"9;(p) is the transformed component of the 
generalized force Q;“ due to the motion in the coordi- 
nate g;; Qi;“(p)q;(p), (7 ¥ 1), are the transformed com- 
ponents of the generalized force Q;” due to the motion 
in all coordinates other than g;; and Q,?(p) is the 
transform of Q;?(t). 

Eq. (27) represents a set of m simultaneous linear 
algebraic equations in the unknown functions 4;(p) . . . 
in(p). 
tule, and the solution may be expressed as 


These equations may be solved by Cramer's 


n 


- P) = 
gi(p) = )(p) = 
jit 2 A(p) Qs?(t 


Qi;( 


A (p) F(p) (28) 
A(p) 

where A(p) is the determinant of the matrix of coeffi- 
cients of G(p) .. . Gn(P); iz is the cofactor of the term 
in the ith row and jth column; and F(p) is the trans- 
lorm of the disturbing force F(t). 

The process of determining the response q;(¢) re- 
quires that the inverse transform of ¢,(p) be taken. 
This may be 
integral? 


accomplished by the convolution 
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(Ap) . ..\ f jAi(p)) 

i(t) = £ I. Fi (= i 1. ‘ 

q( LA(p) P) 0 \A(p) J 
F(t — r)dr (29) 


The inverse transformation of A;,(p)/A(p) can be car- 
ried out by the application of Heaviside’s partial frac- 


tions expansion. This may be written 


JAlb)t _ 4 Addy) 


-{( = 9 PK 
La(p) J Pa A'(p,) 


> 


(30) 
where 
A(p) = (p — pi) (p — fo)... (P — pn) 


and A’(P,) indicates the derivative of A(p) with respect 
to p and the subsequent insertion of the root P,. 

It is seen that the denominator polynomial must be 
factored. This feature of the operational solution is 
one of its most unfortunate aspects, since the factoring 
of a polynomial of high degree with a number of com- 
plex roots is usually an exceedingly difficult and cum- 
bersome task. 

Combining Eqs. (29) and (30) yields the response of 
the 7th generalized coordinate 


. Ais f° . 
— | F(t — r)e’""d 
2d, A’(P,) Jo ' _ 


(2) Solution by Mechanical Admittance 


qi (t) = (31) 


The term ‘‘admittance,’’ as defined in electrical engi- 
neering, is the ratio of the amplitude of the current in a 
circuit to the amplitude of an applied voltage which 
varies sinusoidally with respect to time. In an anal- 
ogous manner, we may define a mechanical admittance 
as being the ratio of the amplitude of a displacement of 
a mechanical system to the amplitude of a sinusoidally 
varying force that causes the displacement. The me- 
chanical admittance is, in general, a function of the fre- 
quency of the applied force. 

In the case of conservative or undamped systems, 
the response is always in phase, or 180° out of phase, 
with the applied force. In the case of systems in which 
damping forces are present, the response is generally out 
of phase with the applied force. In order to relate 
properly the response to the forcing function, the me- 
chanical admittance should therefore be considered to 
consist of two quantities: the ratio of the amplitudes 
and the phase lag of the response with respect to the 
forcing function. 

A comprehensive presentation of the concept of me- 
chanical admittances and its application in the solu- 
tion of vibration problems is given by Duncan.‘ The 
term “‘mechanical admittance’”’ as used here is not re- 
stricted, in its application, to a forcing function and a 
response in the form of a force and a displacement, re- 
spectively, but is considered to apply to any kind of 
response. For example, the forcing function may be a 
sinusoidal force applied at a gun mount, a sinusoidal 
gust velocity pattern, or a sinusoidal control movement; 
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Fic. 2. Unit-step function. 


and the response may be the displacement of a normal 
mode or the stress at a point in an airplane wing. The 
idea for the present extension has been borrowed from a 
development by Dr. C. S. Draper and his associates in 
the Instrumentation Laboratory at the Massachusetts 
Institute of Technology.°® 

Methods for finding the transient response to an 
arbitrary forcing function, from a knowledge of the ad- 
mittance as a function of frequency, have been de- 
veloped and applied in the field of electrical engineering. 
These methods generally involve the evaluation of a 
Fourier integral.2 This is not practical when the ad- 
mittance cannot be expressed as a simple analytical 
function of the frequency, especially since the integra- 
tion must be carried out between infinite limits. A 
convenient approximate method for finding the trans- 
ient response from the admittance function has re- 
cently been developed, however, and is presented be- 
low.® It should be noted that both the exact and ap- 
proximate methods involve the assumption of linearity 
of the system. 

The determination of the mechanical admittance in a 
particular problem may be accomplished conveniently 
by the complex variable method of vibration analysis.° 
If an elastic system is excited by a sinusoidal forcing 
function at a certain frequency, the steady response 
of the system will also be sinusoidal and of the same 
frequency. If the forcing function is represented by 
Fe’ and the response by Ae™, the mechanical ad- 
mittance will be given by the complex number, // = 
A/F. The absolute value of // will give the ratio of the 
amplitude of the response to the amplitude of the forcing 
function, and the argument of //7 will define the phase 
lag of the response with respect to the forcing function. 
that is, we represent the forcing 
me- 


The 


If we put Ff = 1 
function by the unit rotating vector e'“—the 
chanical admittance is given simply by /7 = A. 
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value of // depends in general upon the value of @, $9 
that, if w is regarded as a variable, then // will be , 
function of w. 

The function //(w) can be represented graphically in 
either of two ways. The real and imaginary parts of 
H(w), which are designated R(w) and /(w), respectively. 
may be plotted versus w as shown in Fig. la; or the 
amplitude and phase of //(w), designated // and y, 
may be plotted versus w as shown in Fig. Ib. It cap 
easily be seen that for any real physical system, R and 
HT, must be even functions of w, and J and W must be 
odd functions of w, since the significance of a negative 
frequency is simply that of a vector rotating in the 
opposite sense. 

The transient response to an arbitrary forcing fune- 
tion cannot be determined directly from the mechanical 
admittance function in a convenient manner. Hoy- 
ever, the transient response to a unit-step forcing fune- 
tion can be determined, and this information can be 
used to determine the response to an arbitrary forcing 
function by the application of the principle of super 
position. The unit-step function is shown graphically 
in Fig. 2. 

In the approximate method, the transient response to 
the unit-step forcing function is not found directly, 
Instead, the periodic response to a square-wave forcing 
function of the form shown by the solid line in Fig. 3 is 
determined. If the system is well damped and the 
wave length of the square-wave is taken sufficiently 
long, the response will tend toward an asymptotic 
value, which is the static response, as the end of each 
half-wave is approached. This is shown by the broken 
line in Fig. 3. Under these circumstances, the response 
to the first half-wave for 0 < ¢ < (w/wo) can be regarded 
as al approximation to the response to a unit-step func- 
tion. It can be seen that the response function must 
initially be zero. 

The determination of the response to the square- 
wave forcing function is carried out as follows. The 
square-wave function is represented in the form of a 
Fourier series, 


+ (2/2) > (sin mwot/n) 
n=1,3,5 

that is, it is considered to be a superposition of sin 

soidal functions. The periodic response to the square: 

wave forcing function may then be considered to be a 


superposition of the sinusoidal responses to the sinu- 
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widal forcing functions of which the square-wave is 
Since the first term of Eq. (32) isa constant, 
the response to it will be the static response—that is, 
H/2, where J/\) is the value of |/7; atw=0. The 
response to the function Sin ” wy t can be obtained from 
Fig. 1b and is written |/Z), sin (nwt + y,), where | /7|, 
and y, are the values of |H7| and y, respectively, at w = 
The net response to the square-wave function is 


composed. 


Nw. 


now written 
A(t) = (1/2) | Ho + 


(2/7) = [\/7, sin (nat + y,)/n] 


%=21,3,9... 


(33) 


This equation may be written in the alternative form 
A(t) = (1/2)Ro + 


(2/m) > 


[(R,, sin nwt + I, cos nat)/n}] (34) 
1,3,5 . 


If the system is damped, as it must be if the method 
isto be applicable, the value of |/7| will tend to zero as 
wbecomes large. This fact, together with the presence 
of nin the denominator of the terms in Eq. (33), leads 
to rapid convergence of the series. In consequence of 
this, a satisfactory approximation will usually be ob- 
tained by taking only a moderate number of terms in 
the series. 

When the approximate response to the unit-step 
function has been obtained, the response to an arbi- 
trary forcing function, F(t), is given by the superposi- 
tion theorem? 


t 
J F(r)A’(t — 7) dr 
J0 


where A’(é) is the first derivative of A(é) with respect 
to €. 

The method is now illustrated by application to the 
Solutions of the form ge" are 
into the differential 


x(t) = (35) 


solution of Eqs. (24). 
assumed. These are introduced 
equations to yield a set of algebraic equations of the 
form 


} 9 9 


{(o? — w?) M, — Oi: (w) tdi is 


DL Vis" (w)d; = O:?(w) 


jI#t 


(36) 


The quantities Q;;“, Qi;“, and Q;” are analogous to 
the similarly represented quantities in the Laplace 
transform solution [Eq. (27)], except that they are now 


lunctions of the forcing frequency w. The quantity 


Q,” is determined for a forcing function e'“’—that is, a 
simple harmonic function of unit amplitude. 

The unknowns g;, which are mechanical admittances 
of the system, may be determined by a solution of the 
" simultaneous linear algebraic equations [Eqs. (36)] 
and will themselves be functions of the forcing fre- 
quency w. They will, in general, be complex quanti- 
lies, 

_It is not often possible to determine the quantities 
0, 0., 0,2, as simple analytical functions of w, so 
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that it is not practical to determine the admittance 4; 
in function form. However, if the forcing frequency 
w is given a particular value, the set of equations can be 
solved numerically, and the values of the mechanical 
admittances for that frequency can be obtained. By 
doing this for a series of values of w, the curves of Fig. 
la or Fig. 1b may be obtained by plotting. 

The final step, involving the determination of the 
response of each normal mode to the arbitrary forcing 
function F(t), may be accomplished by the appli- 
cation of Eq. (85). Since, in general, both F(/) and 
the response to the unit-step forcing function will 
probably not be known analytically but only as plotted 
curves, an analytical evaluation of Duhamel’s integral 
will not be possible. This evaluation must, there- 
fore, be carried out numerically or by a graphical proc- 
ess.’ 

It should be noted that the mechanical admittance 
method can often be employed where the Laplace 
transform method would be entirely impractical, either 
because of inversion difficulties or because of difficulties 
in the transformation of the equations. A more 
rigorous mathematical development of the method, 
through the Fourier integral, is presented in Appendix 
(B). 

THE TRANSIENT STRESSES 


(VI) DETERMINATION OF 


If the deformation of an elastic structure has been 
determined, it is always possible, in principle, to deter- 
mine the stress distribution corresponding to this def 
ormation. If the elastic characteristics of the structure 
are available in the form of influence coefficients or in 
fluence functions, the loads required to hold the struc 
ture in the given deformation configuration may be 
computed. Using these loads, a stress analysis of the 
structure may be carried out in the conventional man- 
It is also theoretically possible to determine the 
stresses the with the 
curvature of the elastic structure. 

Although straightforward in principle, these meth 
ods may not be the most direct, particularly if the def 
ormation configuration has been computed in terms of 
Two alternative 


ner. 


from strains associated space 


displacements of normal modes. 
methods for the determination of the transient stresses 
are now considered. 

In the first it is observed that, because of the assump 
tion of small deformations, the principle of superposi- 
tion applies, and the total stress pattern for any def 
ormation of the structure may be considered to be 
a superposition of the stress patterns corresponding 
to displacements of each of the normal modes. There- 
fore, the stress at a particular point 7 in the structure 
may be expressed by a relation of the form 

1 
o, = Lu 


r=m+1 


where the quantity A,” is a constant that represents 
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the stress at the point j due to a unit displacement in the 
rth normal mode. An illustration of the determination 
of A,” in a simple case is given in Appendix (C). 

The stress corresponding to a unit displacement in a 
normal mode is the same as the stress produced by the 
inertia loading associated with a free vibration of unit 
amplitude in that mode. This inertia loading may be 
expressed in terms of the mode shape and frequency of 
the mode, and the constant A;‘” may subsequently be 
determined by a conventional analysis. 

In the second method of approach, the stresses that 
would result if the system were restrained against vi- 
brating are first calculated, and the additional stresses 
resulting from the vibratory motion are added. The 
stresses computed in the first step, assuming no vibra- 
tion, are those that are computed by the airplane de- 
signer when he assumes that the structure is rigid; 
they wili be referred to as static stresses. The addi- 
tional stresses resulting from the vibratory motion 
represent the dynamic overstress that is a consequence 
of the lack of rigidity of the structure. This method 
(which is particularly suited, but not restricted, to un- 
systems subjected to prescribed external 
yas first suggested by Williams.* 


damped 
forces) 

The method is first developed for an undamped elastic 
system under the influence of prescribed external forces. 
In the equations of motion for the vibratory modes 
[Eqs. (21)], the vibratory inertia effect may be elimi- 
nated by setting g, = 0. This yields the static mode 
displacements 

Gr (static) Q,(t) Mw," 

Substituting this expression into Eq. (37), the static 
stress at the point 7 becomes 


> A;” [0,-(t)/M,?] (39) 


F 5 (static) Pepe =i 


where the summation includes only terms associated 
with the vibration modes. It can be shown that the 
stresses given by Eq. (39), if all the normal modes of 
vibration are taken into account, is the same as the 
stresses that are commonly computed by the airplane 
designer when he makes the assumption of structural 
rigidity with the externally applied forces put into 
equilibrium with rigid body inertia forces. 

The total displacement of the rth vibration mode, 
including vibratory effects, obtained from Eq. (21), 
would be written 

gr = [0,(t)/M,@,"] — (G,/w,?) (40) 
or 
dr = Qr(static) (q, ‘w,”) 


Substituting Eq. (40) into Eq. (37) and introducing 
Eq. (39), the total stress at the point 7 becomes 
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— (r) » 
a * tistatic) (A; W, Ur (4] 
r m+1 


This equation shows that the stress in the structure a 
any instant is equal to the static stress at that instay; 
plus an additional stress that is a linear function of the 
accelerations of the vibration modes. In this Case, the 
accelerations of the vibration modes can be derived fron 


Kq. (22) 
t 
Gr, = (/M,) J cos w,(t — 7) (d/dr) Q,(r) dr (49 
0 


It should be noted that the so-called ‘‘static’’ stress jg 
really a pseudo-static stress, since it includes the effect 
of the inertia forces associated with rigid body accelera. 
tion. 

When the possibility of aerodynamic forces or, jn 
general, forces other than elastic or inertia forces that 
depend upon the motion of the system is allowed, the 
solution is derived from Eq. (24), and the stress may 
be written . 


n 
0; _ Ti static) + p is A;(Q,™ M,w,") ie 
r=m+1 
n 


} ® (A, 


r=m+1 


9 , 
@,")Gr 43 


where the first term represents the stress that would be 
obtained under the assumption of rigidity; the second 
term represents the stress due to forces, other than 
inertia forces, which result from the motion of the sys- 
tem; and the last term represents the stress due to the 
vibratory inertia forces. Since each Q,” may be a 
function of all the coordinates, it can be seen that the 
second term may become exceedingly complicated. 
Both methods give equivalent results if enough terms 
in the series are considered. However, for practical 
purposes, only a relatively few terms may be taken 
It is therefore desirable to use, whenever possible, the 
method giving the more rapid convergence of the series. 
In Eq. (41) the components corresponding to the static 
application of the loads are all contained implicitly m 
the first term, while in Eq. (37) they are contained i 
This suggests a more rapid 
In fact, in many 


each term of the series. 
convergence of the series in Eq. (41). 
practical problems involving airplane structures, vibra- 
tions of only the lowest two or three modes are excited 
significantly, and, therefore, if Eq. (41) is employed, 
only two or three terms of the series need be included. 
On the other hand, if Eq. (37) is employed, a sufficient 
number of terms must be taken to ensure that all signil- 
icant components of the static deformation have bee! 
included. 

In spite of the advantages of convergence in thi 
second method, it is unwieldy in its application whet 
Q:™ is not equal to zero—that is, when Eq. (43) musi 
be used. In some cases of this kind,:the use of Eq 
(37) and a larger number of modes may prove to ifl- 
volve less work. 
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AIRPLANE C.G. 


z 


Axes of unrestrained airplane 


Appendix (A) 


FREE VIBRATIONS OF AN UNRESTRAINED ELASTIC AIRPLANE 


The problem of free vibrations of a continuous unrestrained elastic system can be formulated in terms of a 


homogeneous linear integral equation. 
problem of free vibrations. 


cally as the displacement at one point of the system due to a unit load applied at another. 


The kernel is equivalent to the Green’s function of the boundary value 


For a restrained elastic system, the Green's function can be easily interpreted physi- 


For an unrestrained 


system, the Green's function can be defined as the displacement of one point with respect to axes through the 


center of gravity of the deformed system due to a unit load applied at another. 


The unrestrained system is con- 


sidered to be in equilibrium such that the unit external load is balanced by the inertia forces corresponding to the 


acceleration of the system as a rigid body due to the unit load. 
In the case of an unrestrained airplane (Fig. 4), a set of x, y, s axes is chosen to coincide with the principal axes 


of the airplane. 
y directions. 


The airplane is assumed to have small dimensions in the z direction compared to those in the x and 


Because of a unit load in the z direction at (¢, 7), the system undergoes an acceleration composed of the trans- 


lation of the center of gravity and the rotation about the x and y axes. 


written as 


a(r, s) = 


(1/M) + (ér 


The acceleration at any point (7, s) can be 


I,) + (ns/Iz) (A1) 


where .\/ is the total mass of the airplane and J, and J, are its moments of inertia about the x and y axes, respec- 


tively. 
f(r, s) = 


l 
ais 
M 


where p(r, s) is the mass distribution. 


The inertia force due to an element of mass at (7, s) is 


+ *) p(r, s) dr ds (A2) 


I; 


In obtaining the displacement of the system, we find first the displacement with respect to a reference plane that 


is tangent to the elastic curve at the origin. 


Thus, when G(x, y; & 7) is defined as the displacement in the z direc- 


S) 


tion of a point at (x, y) due to a unit force at (£, 7) when the origin of the x, y, z axes is clamped, the displacement 
of a point (x, y) with respect to the reference plane due to a unit force at the point (£, 7) and the distributed inertia 


, : ‘ ] tr ns 
G(x, y; &n) — J f G(x, ¥; 7, S) | + + | p(r, s) dr ds 
‘ i M  P Ie 


force is 
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The integration is taken over the entire surface representing the elastic airplane. If the plane containing the ty) 
principal axes of the deformed system is taken as the reference plane, the displacement at (x, y) is 


K(x, y; & 0) = G(x, y5 & 1) — J J G(r, s; & 9) E + 24+ dl p(r, s) dr ds — 
M i. Fe 
: ] gr ns 
Gls, ¥; #,.8) oa a (ry, s) dr ds + 
Jf ( 3 E l, ™] p ar ¢ + 
JS JS ee J l 4 tp + xr 4% + ys 4 id 4 ™gr + typs 4 m4 . 
VE MI, MIL, 1,? Ll, i? 


p(p, g)p(r, s) dp dqdrds (43 


This is the Green’s function of the unrestrained airplane. The first two terms on the right-hand side of the equa. 
tion represent the displacement with respect to the plane that is tangent to the elastic curve at the origin, and the 
last two terms locate this plane with respect to the plane containing the two principal axes of the deformed sys. 
tem. The function G(x, y; & 7) is symmetrical because of Maxwell's theorem of reciprocal deflections. It can be 
seen from Eq. (A3) that the Green’s function A(x, y; & 7) is also symmetrical. 

The homogeneous linear integral equation that governs the free vibrations of the unrestrained elastic airplane 


2(x, y) = w? JS K(x, y; &, n)p(é, n)2(&, 7) dé dn (A4 


where 2(x, y) is the amplitude of vibration of the oscillating surface. 
According to the Hilbert Schmidt’ theory of integral equations, there are associated with Eq. (A4) an infinite 
number of characteristic values or natural frequencies, w;, and characteristic functions or mode shapes, ¢;(x, y). 
In addition to the infinite number of mode shapes and frequencies associated with Eq. (A4), there are three 
modes of zero frequency. These represent solutions for the airplane as a rigid body. The corresponding fune- 


becomes 


tions ¢;(x, y) can be expressed as follows: 


di(x, y) = C (A5 
g(x,y) = Cox (A6 
o3(x, y) = Czy (A7 


where the C’s are undetermined constants. 
The last two integrals in the expression A(x, y; & 7) [Eq. (A3)| contain either the zero power or the first power 
of andy. These two terms will not contribute anything in Eq. (A4), since the integrals 


/ fv n)2(&, n)dé dn, / / p(é, né2(k, n) dé dy 
J / p(é, n) n2(&, n) dé dy 


are zero. By dropping these two terms, the integral equation of the free vibrations of an elastic airplane may be 


2(x, y) = wf | K'(x, v; & n)p(&, n)2(&, n)dé dy (AS 


K'(x, y; &n) = G(x, Ss 25) = J / G(r, s; &, n) E re : 4 *] o(r, s) dr ds (A9 


In a complex airplane structure, approximate solutions to Eq. (AS) must be devised. The two methods of ap- 
proximation described in Section (II) may be used, both of which lead to a set of linear algebraic equations. In 
the first, the integral equation is replaced by its equivalent set of linear algebraic equations’? 


written as 


where 


iL Dij2j = 21 ‘ws? (A10 
] 
where 


Se Z y . 5 : 
Dy = Cu —- — )Cym: - 7 > Cymy. — my > Cymer; t m (All 
d =e * U7] 
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This is equivalent to dividing the structure into a finite number of discrete masses, and it forms an acceptable 
approximation where the structure has numerous concentrated masses and discontinuities. The subscript 7 refers 
to the ith mass m,. The elastic constants C,,, termed flexibility influence coefficients, represent the values of the 
function G(x, y; & ») computed at the points 7(x, y) and j(&, 7). Therefore, the quantity C;,; represents the dis- 
placement of the 7th concentrated mass in the s direction due to a unit force at the jth concentrated mass when the 
structure is clamped at the origin of the x, y,zaxes. In the second method of approximate solution of Eq. (A8), the 
deformation of the elastic surface can be considered as the superposition of several assumed deformation shapes 


y,(x, y) which satisfy the boundary conditions, so that 


! 


2(x, y) = > vilx, y)B: (A12) 
+=1 


where 6; is a multiplier that determines the amount of contribution of each assumed deformation shape. Substi- 
tuting Eq. (A12) into Eq. (AS) and requiring that the resulting equation be satisfied at a finite number of points 
yields a set of linear equations in 8;._ The characteristic roots and values of this set of equations, as well as the 
set represented by Eq. (A10), can be determined by the process of matrix iteration. 


Appendix (B) A) = 1R(0) ‘i 
2LATION BETWEEN THE INDICIAL ADMITTANCE ANI 1 “7 ) , 
—— ee fg = a ; = = an _ ee Fann sin wt + i) cos wt bdo (B6) 
THE MECHANICAL ADMITTANCE eve t ws o f 


The indicial admittance and the mechanical ad- 


mittance are related, mathematically, by the infinite It is apparent that Eq. (B6) reduces to Eq. (34) if the 


integral is replaced by a s ation. 
integral theorem. !! integral is replaced by a summation 


1 {t° Hw) . ‘ . 

A(t) = = / —e' "dw (B1) Appendix (C) 
Qn 1w 

DETERMINATION OF STRESSES IN TERMS OF 


In general, the integrand has a pole at the origin, and it 
DISPLACEMENTS OF NORMAL MODES 


can be evaluated by residue theory. However, this 
difficulty may be overcome by writing the expression If the deformation of an airplane structure, repre- 
for A(t) in another form: sented by an elastic surface 2(x, y), is known, the load 
| "+© Ti) — H(0) distribution required to hold the structure in the given 

twl ° . . . ° 
A(t) = H(O)1(t) + = / ; eda deformation shape is determined by a solution for 


2a 1a) (Be w(x, y) of the integral equation 
“) 


where the symbol 1(/) represents the unit step function 2(x, y) = / / K'(x, y; & n)w(t, n) didn (C1) 


I(t) = Ofort< O\ 


‘ (B3) 7 > def ats ae heen © ws 2 ——— bo 
1 fort > 0f Do If the deformation has been computed in terms of dis 


placements of normal modes such that 


It can be shown that n 
2(x,v) = >> ¢,(x, ¥) @ (C2) 


M(t) = (1/26) / (e /iw) dw (B4) 


where @,(x, y) satisfies the homogeneous integral equa- 
Transforming Eq. (B2) into real form and putting — tion 


H(iw) = R(w) + iI(w), 


(x, y) = 
l ' F * 
A(t) = 5 RO) + w,? / / K'(x, y; & n) p(& n)o-(& 0) dE dn (C3) 
oe a: I(w) ‘ > de ation shape may be expressed by 
l / fR) ies at Oe one yon (BS) the deformation shape may be expressed by 
2a lw w f 
a(x, y) = 


lf Rlw _ ae ° 42 = oe ¢ n n . > = sea " : 
(w) is an even function [R(—w) R(w)| and I(w — / f K(x, 3. & 0) o-(& n)o(& dé dn (C4) 
r=4 


is odd [J(—w) = —I(w)], Eq. (B5) reduces to 
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Comparing Eqs. (C1) and (C4), there results 
w(x, y) = do o,(x, y)p(x, ¥) w,74, (C5) 
r=4 


Eq. (C5) expresses the load distribution w(x, y) required 
to hold the structure in the given deformation shape 
z(x, y) in terms of the inertia loadings associated with 
free vibrations of the modes. 

As an illustration of the calculation of the stress at a 
particular point in the airplane structure, the shear 
and moment at the section A—B at a distance yp from the 
wing root shown by Fig. 4, are given by the following: 


5; = J / w(x, y) dx dy = 


A 
2, wo) IS $,(x, y) p(x, y)dx dy ; q, = a A ;q, 
A 


(C6) 


where 


a oe J J r(x, ¥)p(x, y) dx dy 
A 
and 


M; = J J w(x, vy) (vy — yo) dx dy = 
A 
Lo, ; / / or(x, y)p(x, ¥) (y — yo)dx dy ; x 
r=4 
A 


= } Bisq, (C7 
r=4 


where B” = off fo,(x, v)p(x, viv — Yo)dx dy 
A : 


In both Eqs. (C6) and (C7) the integrations are taken 
over the shaded area shown by Fig. 4. 
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Comparison of Analytical and High-Speed 
Mechanical-Calculator Solutions of the 
Lateral Equations of Motion of an 
Airplane, Including Product-of- 
Inertia Terms 


M. W. FOSSIER,* R. W. BRATT,7 ano D. G. DILLt 
Douglas Aircraft Company, Inc. 


ABSTRACT 


Two methods are derived for solving the lateral equations of 
motion of an airplane slightly disturbed from steady symmetric 
flight. One method consists of solving the equations analytically, 
and equations are derived for the motion of an airplane resulting 
from constant applied forces and moments. The second method 
gives a step-by-step solution in matrix form, suitable to calcula- 
tion by a high-speed mechanical computer such IBM 
calculator. Solutions computed by IBM calculators may be 
done with a considerable saving in time as compared with the 
This saving is a function of the number of 


as an 


analytical solutions. 
cases to be calculated. In addition, the possibility of error, high 
in the analytical solution, is reduced considerably. 

hypothetical sweptback-wing 


It is found that there is 


Calculations are made for a 
fighter-type airplane by both methods. 
no difference in the results as calculated by each method. 

The effect of including product-of-inertia terms in the calcula 
tions of response characteristics and dynamic stability is investi- 
It is found that, for the longitudinal 
the flight path, the 
stability and 


gated on the same airplane. 


principal axis inclined above inclusion of 


product-of-inertia terms improves the decreases 
the rolling performance from that calculated without the product- 


of-inertia terms 


NOTATION 


o = angle of bank, rad. 

y = angle of azimuth, rad 

8 = angle of sideslip, rad. 

¥ = angle of flight path to horizontal, rad. 

v7] = angle of attack of principal longitudinal axis of 
airplane, rad. 

r) = mass density of air, slugs per cu.ft. 

V = speed of airplane, ft. per sec. 

q = dynamic pressure, lbs. per sq.ft. 

S = area of wing, sq.ft. 

b = span of airplane, ft. 

m = mass of airplane, slugs 

Mm = m/pSb 

Kz, = radius of gyration in roll about principal longi- 
tudinal axis, ft. 

K,, = radius of gyration in yaw about principal normal 
axis, ft. 

I, = mK,,?/qbS 

I,, = mK,,?/qbS 


Received July 6, 1949. 

* Aerodynamicist. 

+ Aerodynamics Engineer 
t Research Engineer. 


I, = ],, cos? + I,, sin? 7 

' = I,, cos? 7 + Iz, sin? n 

Tes = (J,, — I,,) sin 7 cos 7 

CL = lift coefficient (assumed constant) 
Cy = lateral-force coefficient 

Ci = rolling-moment coefficient 
Cr = yawing-moment coefficient 
Cig = 0C,/d8 

Cnrg = 0C,/08 

Cys = oCy/ds 

Ci, = 0C;/0(dp/dt) 

Cng = OC,/0(d¢/dt) 

Cy, = WCy/d(dp/dt) 

Cy = 0C,/d(dy/dt) 

¢ ny = OC,/O0(dy/dt) 

Cy, = dCy/a(dy/dt) 

D = d/dt 

p = Dé 

pb/2V = wing-tip helix angle 

Cip = 0C; O( pb 2V) = (2V b)Cig 


INTRODUCTION 


— RESEARCH'? HAS SHOWN the need for 
including product-of-inertia terms in the lateral 
equations of motion to determine the lateral stability 
of an airplane. It was felt desirable to include these 
terms in the calculated response of a typical sweptback- 
wing fighter-type airplane to constant applied rolling 
and yawing moments and side force. It is the purpose 
of this paper to present the results of these calculations, 
along with two methods of calculation. 

While the analytical procedure is straightforward, 
it requires a considerable amount of hand computation. 
(punched card) 

It was found 


The possibility of using automatic 

calculators was therefore investigated. 
that, as often happens, the methods suited to hand 
calculations were not for the machines. 
The method developed herein employs a direct step- 
by-step integration of the original lateral equations of 


optimum 


motion. 

It should be pointed out that, although the solutions 
apply strictly to infinitesimal disturbances only, they 
have been found to apply with fair accuracy to finite 
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disturbances that are sufficiently large to make the 
conclusions drawn of practical interest. 
EQUATIONS OF 


SOLUTION OF THE 


MOTION 


(I) ANALYTICAL 


(A). Aileron Roll 
The lateral equations of motion for an airplane slightly 
disturbed from steady symmetric flight may be written! 


(I,D* — Ci,D)o + (—122D* — CyD)b — CiB = Cr 
(1) 


[,D? — C,,D 
—I,:D* — C,,D 


—C, — Cy,D 


(2 - c.) 


If the forcing functions are constants, the particular 
solution for ¢@ will be of the form 


04) = Ky't + K,’ (5) 


Thus, the complete solution for ¢ is 


nN 
@= > Ce“ + Ki't + K,’ (6) 
i 1 
where 
tan y¥(—C,C,, + Cigna) 
K,' = : : 


—CiyOn, + CigCny + tan ¥(Ci,Cr, — CigCrg) 


The \; are the roots of the fifth-power polynomial in 
D obtained from Eq. (4). The C; are determined from 
the boundary conditions, and AK,’ may be determined 
from Eq. (6) after the C; are known. It can be seen 
by inspection of A that one of the roots of this equa- 
tion will be zero, thereby reducing the equation to a 
quartic. 

For modern airplanes, there are usually two real 
roots (Ai, Ae) and a pair of complex roots (@ + 7) of 
If \y, Ae, and a are negative, 


the quartic equation in D. 
“dutch- 


the aircraft is stable. Positive a 
roll” instability, and a positive value of one of the 
The quartic 


indicates 
real roots indicates spiral instability. 
from which the roots are obtained may be written 
as 


Ad! + BB+ Ch? + EX + F = 0 (7) 


where 


AERONAUTICAL 


—I,:D° — CyyD 
1,.D? — C,,D 
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(—J,,D* aie Cy, ,D)o > (I,D? Pa, C,,D)v = C,,8 “= Cc 


(—C, — Cy .D)o + } [(2bpu 'V) _— Cy ,|D a 

7 ° t f¢ r ’ ’ 

C, tan yy ¥ + [(2bu/V)D — Cy,|8 = Cy, 63 
Since the operator D may be regarded as an algebrajc 
quantity, the equations may be considered as algebraic 
equations for ¢, y, and 8, and the complementary 
function may be obtained by setting the forcing 
functions equal to zero and expanding the determ; 


nant 


-Ci, 
anit, 
8 == 0 
; 2bu 
D — C, tan y Y i? — Cy. 
A= (2bu VI 
B= —(2bp/V)Ks — Cy,l 
’ 2bu ... : _— . oe . 7 
C = "a (K3 + K;4) + Cy ,Ks il Cy Ks si Cy Ks 
: 2bu _. es , - + 
en 2 Ky — C.Ks — C, tan yK3 — Cy, Ky 
Cy A; — Cy, Ke 
F = C,K; — C, tan yKe 


in which 
, = EY, — i 
Ci, + Cn, 


Ky = 1,C,,+ LC, 
Ks = InCi, + I:Ca, 
Ks = C1,Cx, — C..Cy, 
K, = IC, + 11, 
Ke = Cu,Ci, — Ci.Cu, 
K; = Ci,Cay — Cr,Cr, 
Ks = IK: + Ks 


The roots (Ai, Az, @ = in) may be extracted easily by 
synthetic division. 
Using these roots, the complete solution becomes 


@ = Ce + Coe + e (C3 cos nt + Cysin nt) + 
Ky't + ks 'e) 


The coefficients C; of Eq. (7) are determined from the 


boundary conditions at ¢ = 0. The conditions are 
¢ = (Dé) = (Dé).) 
y = ¥, (Dy) = (Dy)o7 4) 
B = Bo \ 





wl 


in 
of 


(D 
(D 


By 
Eq 
anc 


wel 
sist 
ma) 


ebraic 


ebraic 
ntary 
orcing 
termi- 


y by 








SOLUTIONS OF 


LATERAL 


EQUATIONS OF MOTION 273 


Consideration of Eqs. (9) in Eqs. (1), (2), and (3) yields the following conditions at t = 0: 


(DB)o 
(D°o), = CU I)[Ky + Ki(Do), + Ku(Dy), + K38,] 
(Dy), = | T)[Kiz + Kis(D¢), + Ku(Dy), + K38.] 
(D*8), 
(D*o), = ( [)[Ku(D*y), + Kyo(D*¢), + K;(DB),| 
(Dy), = CU T)[Ki3(D*), + Ku(D*p), + Ks(DB),] 


ll 


= [1/(2bu/V)NCy, + Crd. + Cr,(D6)o — [(2bu/V) — Cyj](DW)o + Ci. tan y + Cy,8.} —) 


[1/(2bp V)}iCy,(DB). — [(2bu/V) — Cy] (D*y), + C, tan y(Dy), + Cr(Dd), + Cy ,(D*¢) 0} 


f (10) 


(D38), = [1/(2bu/V)]{C.(D*¢), + Cy ,(D*9)o — [(2bu/V) — Cy) (Do + C, tan y(D*y), + Cy,(D*8) 0} 


(Do), = (1 I) [Kio(D*¢), + Ku(D*y), + K5(D?8),] 
(Diy), = (1/D[Kis(D*¢), + Ku(D*y), + K3(D*8),] 


where, in addition to the A’s listed under Eq. (7), 


Re = fy + ity 
Ky = 1,C,, + InCn, 
Ku TCuy + LesCuy 
Ky» = T12Cy + TC no 
Kis = IC ss Tug 
Kis = T,:Cy + ICuy 


dl 


By using the following substitutions: 


MC) = di, oC = de 
aC + nC, = @3, aC, _ C3 = a 


in the first four derivatives of Eq. (8), the following set 
of simultaneous equations are obtained: 


(Do), — Ki’ = a1 + a2 + 4; 


(D*o), = Nid; + Ade + ad3 + Ay 
(D9), = \j’a, + doe + (a? — n*)as + (11) 
Zand, 
(D‘o), = 8a; + Aa + (a*® — 3an?) X 
a3 + (8a — n*)a4 
By substituting Eqs. (9) and (10) into Eqs. (11), 


Eqs. (11) may be solved for the a,, from which the C, 
and Ky are then determined. 
It is seen by differentiating Eq. (8) that 


b = Do = aye™ + ace + e*(a3 cos nt + 
a, sin nt) + Ky’ (12) 


B. Sideslip 
As mentioned above, the roots of the quartic [Eq. 
(7)| apply for each of the motions, so that the solutions 
differ only in the constants of integration, determined 
by the boundary conditions, and the particular integral. 
The particular solution for 6 is of the form 


B= K;' = 
(Ci, om —_ Cny C,,) a tan ¥(— Ci Cno + Cn Cw) 


rr = ; —— (13) 
(=CijCng + CuyCig) + tan (CiCn, — CagCi,) 


and, therefore, the complete solution is 


B = bye + bee + e%(b; cos nt + by sin nt) + Ks’ 
(14) 


The boundary conditions at t = 0 have been found 
in Eqs. (10) and are substituted into the following 
simultaneous equations: 


B, — Ks’ = b + bs + bs 

(DB), = didi + Avbo + abs + hy 

(D2B), = Ar2b, + Ae + (a? — n*)b3 + Zanhy (15) 
(D*B), = )y3d1 + Aokbe + (a* — 3an?)bs + 


(3a°q — n*)d, 


The 5; thus obtained are then substituted into Eq. 
(14) to give the sideslip solution. 


(II) Matrix SOLUTION 


The method developed herein is that of a step-by-step solution of the lateral equations of motion especially 


well suited to calculation by a high-speed mechanical computer such as an IBM Calculator. 


The method con- 


sists of approximating the solution in any one step by a Maclaurin series, using as many terms of the series as 


may be necessary to obtain any desired degree of accuracy. 
The equations of motion may be written in the following form:* 


Kes | Ce Cy C : 
|! i a | 7 . . p | ~ cb ~ bs | ~ 
—I 5 I, - y Cng Cny y ng Cc no 





2b 
| -c, ( “3 iste Cx.) D H + [-—C, — C, tan y — Cy] 


If to the above set of equations is added the identity 


¢ | 


2bu E a 
A "DB = Cy, (17) 
B J 
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? Do in 
>| ° - 4 = 0 (18) (I 


Eqs. (16), (17), and (18) may be combined into the following matrix equation: 






























































Tro 0 o o|lfe] f ® 0 6 <1 0 P A yie 
0 1 0 0 0 y 0 0 0 0 —] y 0 
2bu : ” 5 2bu . 
0 0 Y 0 0 DB 4+]/—-C, —C,tany —Cy, —Cy, y = Cy, B =1C,, 
0 0 O I, —Izs Do 0 0 tr, —Ci, —Ci, Do C, Su 
000 -In I} |D¥} | 0 0 —Cy, —Cy —Cr, Ly) Lc, 
or 
etc 
Te | [0 0 0 I 0 To 7 1) ea 
y 0 0 0 0 1 y obt 
V V V V V 8 
- - - Cy =~ iy Cy, — 1 
2bu Cr 2bu eaiantl 2bu Cr, Qbu *¢ (5 "y ) 
Do 0 0 I (Cif. + Cyl rz) T (C igle + Cn jlzz) I (C ile + Cnylzz) Do 
| h. ‘. = ; 
Dy | 0 0 7 (Ciglse + Cngls) 5 (Cigles + Cugls) 7 (Cyfles + Cagle) || Dv 
™ = - & ai or, 
. 0 7 | 
0 Qn+ 
Vv 
2bu i 
| (19) 
] (Cy, I, ++ C, fs wh 
1 ‘ 
] (Cotes + Cats) ver 
r 7 out 
Eq. (19) may be written symbolically as ; : 
00 
Dq = Aq+ B (20) r00 
is ¢ 
where a 
T rn) of | 
yp nut 
q=| B . 
Do cal 
Dy Eq. 
= | inte 
0 0 0 1 : (29 
0 0 0 0 I 
J j V J stal 
aa - C, te - Cy = Cy, 
a 2bu ° 2bu sieeadl, 2bu ¥8 2Qhu (6, tim 
7 1 1 I ma: 
Ks K 0 (K 
0 0 ' iis peo 7 Aw cak 
. K (K43) K . 
0 0 I \3 7 “" 7 (Kua) ' 
0 whi 
0 T 
B =| (V/2bu)Cy, 
(1/T)(Ko) 
(1/I)(K 











(18) 


(19) 


(20) 











SOLUTIONS OF LATERAL 


in which the K’s are the same as those given in Section 


(I) and 
I = 1,1, — I,,* 
Expansion of g in a Maclaurin series about ¢ = 0 
yields 
At ( At)? ( At)3 


g= G+ 5 dt 5 Gt 


a Got... (21) 
oO. 


Successive differentiation of Eq. (20) gives 


A%g + AB 
Atg + A2B 


q = Ad 
q = Ag 


etc. 


Substituting these quantities into Eq. (21), one 


obtains 
pw [ict cana 4 ary are Jy 
| (ane - = A+ oe a? + ... |B 
or, expanding in a Taylor series for the (m + 1)st step, 
deri = E + (At)A + a ae + ou A+ Ja + 
| (an + a , — At + .. | B_ (22) 


where J; is the unit matrix. 

The time interval would usually be chosen for con- 
venience in plotting and the series expansion carried 
out far enough to give the accuracy desired. However, 
it should be noted that At X Amar. should not exceed 
3 or 4, where Amar. is the largest expected characteristic 
root. The reason for this is simple, since the process 
is equivalent to the expansion of several exponentials 
in a series, and it is well known that the direct expansion 
of e* for large values of x requires an extremely large 
number of terms. 

Thus, knowing the conditions at ¢ = 0, one may 
calculate g as a function of time by application of 
Eq. (22) to successive steps taken at any desired time 
interval. Calculations of the type involved in Eq. 
(22) are relatively simple for IBM machines. 

It may be pointed out that the coefficients of the 
stability quartic [Eq. (7)], and thus the period and 
time to damp to one-half amplitude of the motion, 
may be obtained as a by-product of the time-history 
calculations. Expansion of the characteristic equation 
M — A = 0 yields the equation 


A‘ + DAF + A? + dA +e = 0 
which is just Eq. (7). 
The coefficients }, c, d, and e are given by 


b = —1(S)) 
¢ = —(1/2)(S2 + bS;) 
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TABLE 1 
Summary of Stability Parameters for Hypothetical Sweptback- 
Wing Fighter-Type Airplane 


Flaps Up 
Altitude = 10,000 ft. 


CL 0.3 0.6 0.9 
2bu/V 4.17 5.90 7.24 
Cy — 0.00826 —0.01106 —0.01316 
Cy 0.00528 0.00881 0.01273 
Cg 0.002560 0.001103 —0.001447 
Cry — (0.01403 —0.01838 — 0.02048 
Cr, 0 0 0 

Cry 0 0 0 

Cig —0.1454 —0.1838 —0.2248 
Crg 0.2178 0.2077 0.2057 
Crg —0.762 — 0.767 —0.799 

| # 0.003875 0.00782 0.01323 
a 0.03660 0.0730 0.1080 
 - 0.001040 —0.00315 —0.01278 


d = (1/3)(S3 + dS, + cS1) 
—(1/4)(S; + BS; + cS, + dS;) 


— 
where 
S, = > diagonal terms in matrix A 
S; = >> diagonal terms in matrix A? 
S; = >> diagonal terms in matrix A® 
S; = > diagonal terms in matrix A‘ 


This scheme for the determination of the coefficients 
of the stability quartic is generally known as Leverrier’s 
method and is applicable to systems of equations of any 
order. A check should be made of the number of 
significant figures retained in the last coefficients of the 
stability equation, because at times they are obtained 
as the differences of large numbers and all the significant 
figures are lost. If such is the case, the last terms 
may usually be found from evaluating Eq. (7). 

It may be pointed out that the general method may 
also be applied to solve the longitudinal equations of 
motion, as well as systems of equations of any number 
of degrees of freedom with constant coefficients. 


APPLICATION OF METHODS TO SWEPTBACK-WING 
FIGHTER-TYPE AIRPLANE 


(IIT) 


For the purpose of illustrating the methods just 
described, calculations have been made on a hypo- 
thetical conventional fighter-type airplane with a 
sweptback wing. 





























C.*09 AT 10,000 FEET | 
Ixz *0] | | | 
40 2 ! ) ¢—1e 
n ANALYTICAL SOLUTION . 
(DEG/SEC,) —_ 
20 + ——+}+ a 
- : | 
2 LA 4 jek ee 
DEG. | | oon 
re) 1 by aad 
0 5 io 1 20 #25 #30 35 


TIME , SECONDS 


Fic. 1. Comparison of analytical and step-by-step solutions. 
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Fic. 2. Effect of the product of inertia on damping of the 


lateral-directional oscillation of a typical sweptback-wing fighter- 
type airplane. 


10,000-ft. 
altitude were estimated by the methods of reference 4 
Fig. 1 presents a com- 


Values of the stability derivatives at 


and are presented in Table 1. 
parison of the time history of airplane motion as 
calculated by the analytical method and by the ap- 
proximate method using IBM calculators. The agree- 
ment between the two methods is within plotting 
accuracy. 

Calculations have been made to show the effect of 
the product-of-inertia stability 
and response characteristics of the hypothetical air- 
plane. It can be shown that the period and time to 
damp to one-half amplitude of the lateral oscillations 


terms on dynamic 


are given by 
P = 
T, _ = 


2r/n 
—(0.693/ a 


where a and 7 are the same as described in Section (I). 
It is seen from Fig. 2 that the apparent dutch-roll 
instability at high angles of attack indicated by 
neglecting product-of-inertia terms actually does not 
exist. The effect of the product-of-inertia terms is 
pronounced at high angles of attack because the 
longitudinal principal axis of the airplane is inclined 
at a large angle to the flight path. It is also seen that, 
for the longitudinal principal axis inclined below the 
flight path, the stability actually is less than would 
be indicated by calculations neglecting product-of- 
inertia terms. 
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Ixz TERMS INCLUDED IN CALCULATIONS 
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ALTITUDE = 10,000 FEET 
160 Se 


Vj 1:03 





120 fF 


ANGLE OF BANK 
$ , (DEG) 


80 


40 








60 





40 


ROLLING VELOCITY 
fe , (066/SEC) 








ANGLE OF 
SIDESLIP 
A. (0€G) 





fe) 
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Fic. 3. Effect of product of inertia on time history of motion 
of typical sweptback-wing fighter-type airplane; motion pro 
duced by an abrupt 10° (total) aileron deflection with rudder 
fixed. 


This same effect is seen in the time histories of the 


response to an applied 10° aileron deflection and a 


fixed rudder, shown in Fig. 3. It is seen that, at angles 
of attack in which the longitudinal principal axis 1s 
inclined above the flight path, omission of the product- 


(Continued on page 296) 
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Some Recent Measurements in a JT wo- 
Dimensional Turbulent Channel 


JOHN LAUFER* 
National Bureau of Standards 


SUMMARY 


The paper describes detailed measurements in mean and turbu- 
lent velocity fields in a two-dimensional fully developed turbulent 
channel flow, including distributions of velocity fluctuations, 
correlation coeffictents, and microscales. It presents results of 
measurements in the laminar sublayer and discusses Reynolds 
Number effects on the turbulent velocity field. Measured 
spectra of the velocity fluctuations and of the Reynolds shear 
permit some conclusions regarding the locally isotropic structure 


‘ 


of the channel flow. 


SYMBOLS 


x = longitudinal coordinate in the direction of 
the flow (x = 0 corresponds to the chan- 
nel exit) 

y = lateral coordinate (y = O corresponds to 


the channel wall) 


lateral coordinate 


a Oe 4 = distances (in the x, y, 2 direction, respec- 
tively) between points at which wu’ cor- 
relations are measured 

d = half width of the channel (2.5 in.) 

u = mean velocity at any point in the channel 

Ur = maximum value of the mean velocity 

u'= V x” = velocity fluctuations in the direction of the 
mean flow—i.e., x 

y = V 02, w’ V w’? =velocity fluctuations normal to the 
mean flow in the direction y and z 

u'/u V u’2/u2, v'/u V v2 u?, w'/u V w’ u2 

p = pressure at any point in the channel 

p = air density 

T = total shearing stress 

T0 = shearing stress at the wall 

Ur = V t0/p 

M = absolute viscosity of air 

! = kinematic viscosity of air 

R = Reynolds Number based on half width of 


channel and maximum mean velocity 


correlation coefficient responsible for the 
apparent shear 


microscales of turbulence; 1/A,2 = 


oe) 2 me) 
(= wae hy? S (= aig 
2 ok) 
To (F _ 


Based on a paper presented at the Fluid Mechanics Session, 
Seventeenth Annual Meeting, I.A.S., New York, January 24-27, 
1949. Received November, 1949. 

* Physicist. Formerly Research Fellow, Guggenheim Aero- 
nautical Laboratory, California Institute of Technology. 


u’'(O) u’(X) 
Rr = , Ry = 
u’(O) u’(X) 


u'(O) u’(Y) 
u’(0) u’(Y)’ 
_ u’(0) u’(Z) 


R = r, 
u'(O)u'(Z) 


n = frequency (cycles per sec.) 


F(n), Fy2(n) = fraction of the turbulent energy u’ associ- 


ated with a band width dn (normalized 


form) 
Fy'p’(n) = fraction of the turbulent energy shear 
u'v' associated with a band width dn 


(normalized form) 


INTRODUCTION 


— PURPOSE OF THIS INVESTIGATION is to obtain ex- 
perimental information concerning the mechanism 
of nonisotropic turbulence. Because of the compli- 
cated nature of this phenomenon, it is essential to carry 
out the investigation in flows having the simplest con- 
figurations—that is, in flows that can easily be treated 
theoretically and at the same time are amenable to de- 
tailed measurements. 


The simplest type of shear flow of course is the 
Couette flow. However ;this requires a moving bound- 
ary, which is difficult to realize experimentally. Pai 
reported measurements in the turbulent flow field be- 
tween two rotating cylinders.’ He has shown that un- 
der such conditions a truly two-dimensional turbulent 
Couette flow cannot be established experimentally be- 
cause of the presence of secondary motions. 

The problem of the turbulent boundary layer is more 
difficult to treat theoretically because of the nonintegra- 
bility of the equations of motion. The difficulties 
encountered in the experimental phase of problems are 
well described in a paper by Schubauer and Klebanoff.” 

The turbulent wake behind a cylinder represents 
another configuration where a nonisotropic turbulent 
flow can be established. Some recent measurements of 
Townsend,* however, show that, except close to the 
core of the wake, the flow is turbulent only for a fraction 
of the total time of observation; the flow is then not 
Similar observations were made by Corrsin in 
Furthermore, the fact that both 


steady. 
a round free jet.* 
mean and fluctuating quantities are decaying with dis- 
tance presents additional complications. 

In the case of the turbulent channel flow the equa- 
tions of motion are simple, and the boundary conditions 
defined. Earlier investigators®*® apparently 


are well 
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encountered no serious difficulties in the experimental 
realization of such flow. It was therefore decided to 
choose this type of flow configuration for the study of 
the nonisotropic turbulent problem. 

The present investigation was carried out under the 
sponsorship of the National Advisory Committee for 
Aeronautics at the Guggenheim Aeronautical Labora- 
tory, California Institute of Technology, and is reported 
in more detailed form in reference 7. The author 
wishes to acknowledge the valuable advice of Dr. H. W. 
Liepmann throughout the investigation. 


TURBULENT CHANNEL FLOW 


The important question that had to be cleared at the 
beginning of the investigation was: Can the theoretical 
assumpttons of fully developed turbulent flow be re. 
alized experimentally? This question is not trivial as 
is often believed, because there exists a possibility of 
secondary flow in the channel and a possibility that the 
fully developed region may not be reached without ex- 
cessive length of channel when the channel must be wide 
enough to permit detailed measurements of turbulence 
- level, scale, etc. 

Indeed, difficulties did occur during the investiga- 
tion. At first a channel of 1-in. width was constructed 
in order to ensure two-dimensional flow (aspect ratio = 
60:1) and to have a large entrance-length width ratio 
(78:1), thus keeping the channel length still within rea- 
sonable dimensions. Measurements have indicated 
that no secondary flow effects were present and that a 
fully developed turbulent region existed at the section 
where the measurements were carried out. However, 


the microscale of turbulence was so small in the channel 
that the length effect of the hot-wire* required large 
corrections to the measurements, and thus the accuracy 
of the results was considerably reduced. It was neces- 
sary therefore to increase the channel width. The final 
measurements were carried out in a 5-in. wide, 23-ft. 
long channel, and it was necessary, of course, again to 
check the two-dimensionality and fully developed na- 
ture of the flow. 

As mentioned above, the simple geometry of the 
channel flow allows an integration of the Reynolds 


equations. The equations have the form® 
Op o (5°) 
y— =7-+ constant; = 0 (] 
~ Ox Ov \Ox 


In the center of the channel the shear vanishes; hence 
if the channel has a width 2d, we have 





d(0p/O0x) = constant 
therefore 
Op du = ; 
y—d =r =pyu— — pu’y’ (2) 
( } Ox dy , 


The shear thus can be obtained from the pressure gradi- 
ent, 0p/dx; since the variation of 7 is linear across the 
channel and u’v’ = 0 at the wall, from the mean veloc- 
ity gradient at the wall, 7) = u(du/dy); or, finally, from 
direct measurements of u’v’, since except very near the 
wall u(du/dy) < pu'v’. 

* The length effect of a hot-wire corresponds to the fact that 
the hot-wire averages fluctuations over its length and that the 


results therefore require correction. 
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MEASUREMENTS AND RESULTS 


(1) Reynolds Number Range 


The investigation was conducted in a channel of 12: 1 
aspect ratio which was attached to the two-dimensional 
tunnel at GALCIT as shown in Fig. 1. The detailed de- 
scription of the tunnel, the channel construction, and the 
dimensions are given in reference 7. The experiments 
were carried out at speeds Up = 3, 7.5, and 15 m. per 
sec., corresponding to Reynolds Numbers 1.23 X 10‘, 
3.08 X 104, and 6.16 X 10+. 


2) Pressure Distribution 


Pressure measurements along the channel center are 
given in Fig. 2 for the two higher Reynolds Number 
For the lowest Reynolds Number flow, the 
gradient was extremely small (approximately 0.0003 
mm. of alcohol per cm.), and accurate measurements 


cases. 


were not possible. The figure also gives the pressure 
gradient obtained from the mean velocity gradient at 
the wall 


From Eq. (2)—namely, for y = 0 


—d(Op/Ox) = to = u(du/dy), =o 
Itis seen that the results of the two independent meas- 
urements are in good agreement. 
3) Mean Velocity Distribution 


Fig. 3 gives the mean velocity distribution across the 
channel at x ~ —2 in. for the three Reynolds Number 
tases. Careful measurements were made close to the 


2 4 16 18 


given in Fig. 6. 
crease with increasing Reynolds Numbers as already 
obtained by 


son. 
fluctuations relative to the local mean speeds (u’/x) 
close to the wall, even though the absolute values here 


20 22 2¢é 


Pressure distribution along the channel. 


wall, and Fig. 4 indicates the linear behavior of the mean 
velocity within the so-called laminar sublayer—i.e., the 
layer where u(du/dy) 2 pu’v’. For R = 12,300, the 
measured velocities seem to be too low close to the wall. 
This is attributed to the fact that in this region the 
velocity fluctuations are extremely large (u’/u > 0.30), 
and it may be shown qualitatively’ that the mean veloc- 
ity indicated by the hot-wire is smaller than the true 
mean if the local velocity fluctuations are large; how- 
ever, the behavior of the hot-wire in regions of large 
fluctuations is quantitatively not well known, and no 
corrections were therefore applied to the measurements. 
Fig. 5 shows that the mean velocities follow the well- 
known logarithmic law with the exception of the region 
near the wall for values yU,/vy <30. This critical value 
of yU,/v is in good agreement with previous results. 
The lower left curve in Fig. 5 corresponds to a linear 
velocity distribution, and the approach toward this 
velocity distribution near the wall is clearly seen. 


(4) Intensity of Turbulent Velocity Fluctuations 


The measured longitudinal velocity fluctuations are 
It is seen that the fluctuations de- 
Wattendorf.6 The measurements of 
Reichardt are also indicated in the figure for compari- 
It is of interest to examine the behavior of the 


should be accepted with reserve because the linearized 








SCIENCES 


MAY, 


1950 



















































































280 JOURNAL OF THE AERONAUTICAL 
© =/2,700 
© = 30,800 
© -6/,600 
— ee 
3 ee —— —_ a —. 
Zu | 
iA 
| 
Fe) - te 
aH 
Us 
4 “es -— 
ae — sai — a cn — —— 
O 
OQ es ihe pg 4 ei Ss e 3 3 J 10 
= & 
d 
Fic. 3. Mean velocity distribution. 
| _[. 
—_—" |_-9 
————— | ee 
| "a a J 
Al 4 dé 


Sls 





© R*=/2,300 


® A=30,800 
© 4-=6/,600 


























Mean Veétoci7ry Disreisution Nere THE Wate 












































.O/ 


02 


.O3 


R 
Ss 
a 


06 


O7 


-O9 


JO 


S/ 





2 


le 


|< 











TWO-DIMENSIONAL TURBULENT CHANNEL 281 





LOGARITHMIC REPRESENTATION OF THE MEAN VELOCITY DISTRIBUTION 


© R=!2,200 (1" CHanner) 
© R=12,300 (5" CHannet) 
© R= 30,800 (5* Cnawnec ) 
@ R=61,600 (5° CHANWEL) 
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hot-wire theory is not strictly applicable to large fluctu- 
ations. Fig. 7 indicates that the values of u’/u reach a 
maximum well within the laminar sublayer, the dis- 
tance from the wall depending on the Reynolds Num- 
ber, and then approach a finite value (u’/u ~ 0.18) at 
the wall independent of the Reynolds Number. It 
should be noted that Taylor already pointed out, based 
on the ultramicroscope measurements of Fage and 
Townend, that at the wall u’/u and w’/u approach a 
finite value.® 

The distribution of v’ and w’ is shown in Figs. 8 and 9. 
The X-type hot-wire anemometer was used for the 
measurements of these quantities, and, because of the 
comparatively large size of the meter, no measurements 
could be made near the wall. The values of v’ and w’ 
are nearly the same, w’ increasing somewhat faster 
toward the wall as also observed by Fage and Townend. 


(5) Turbulent Shear Measurements 


' was meas- 


The correlation coefficient k = u’v’/u'v 
ured (Fig. 10) and was found to decrease with increasing 
Reynolds Number. However, comparison with meas- 
urements made by Wattendorf, Reichardt, and the 
present author in the 1l-in. channel indicates that k is 
not only a function of the Reynolds Number but de- 
pends on a characteristic length parameter, such as the 
channel width. 


= 
e 


From the values of u’/ Uo, v’/ Uo, and k the turbulent 
shear can be computed and can be compared to the 
value calculated from pressure or mean velocity gradi- 
ent measurements. In the case of the two lower Rey- 
nolds Numbers the results are self-consistent; for R = 
61,600, however, the turbulent shear obtained from the 
hot-wire measurements is about 20 per cent too low. 
The reason for this discrepancy is not understood as yet. 


(6) Measurements of Spectrum 


The spectrum distribution F(m) of the w’ fluctuations 
was obtained at different values of y/d. A Hewlett- 
Packard wave analyzer was used for these measure: 
ments. In order to reduce signal-noise ratio in the 
high-frequency range, a hot-wire of 0.00005 in. diameter 
was chosen having a time constant of approximately 
0.20 millisee. With this method it was possible to de- 
tect energy values in the high-frequency bands that 
were 10~’ times the value corresponding to that meas 
ured at the lowest frequencies. Because of the limita: 
tion of the wave analyzer, the spectrum below 30 cycles 
could not be obtained. A typical spectrum distribu- 
tion is shown in Fig. 11. It is seen that, in the high- 
frequency region where viscosity plays a dominant 
role, the distribution follows over a wide range the 
law predicted by Heisenberg.’® It should be met- 
tioned, however, that in the region in question the ratio 
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of noise to signal is high; it increases from 10 to 8O per 
cent in the range 2,000 to 5,000 cycles. This fact and 
possible wire length corrections make the results still 
somewhat uncertain. On the other hand, ‘since the 
measurements were well reproducible whether using 
mean square or root mean square output meters, it was 
felt to be reasonable to make comparisons with the 
theoretical prediction. 

By the introduction of the concept of local isotropy, 
Kolmogoroff helped to clarify somewhat the little un- 
derstood structure of shear flows. Speaking in rather 
loose terms, this hypothesis requires that the smaller 
eddies in turbulent flow approach isotropy. It follows 
that beyond sufficiently high frequencies no correlation 
exists between the components of the velocity fluctua- 
It is of considerable interest to verify experi- 
Townsend in a 


tions. 
mentally Kolmogoroff’s hypothesis. 
turbulent wake has shown the existence of local iso- 
tropy.'' He measured essentially the terms (Ou’/Ox)?, 
(1/2)(dv’/Ox)2, (1/2)(Ow’/Ox)? and found them to be 
This method, however, is not a sensitive test 
In addition, for fairly low Reynolds 


equal. 
for local isotropy. 
Numbers as in the present case, a large part of the 
frequencies contributing to the term (Ou’/Ox)*? can be 
shown to be in the nonisotropic region of the wu’ spec- 
trum. Recently, Corrsin used a more direct method 


for testing local isotropy in a free turbulent jet.'* He 
obtained the spectrum of the correlation coefficient k 
and has shown that with increasing frequency the corre- 
lation spectrum decreases and approaches zero faster 
than the wu’ spectrum. 

Essentially the same method was used in the present 
investigation. Instead of the spectrum of the correla- 
tion coefficient, however, the somewhat easier obtaina- 
ble spectrum of the Reynolds stress, u’v’ was measured. 
Fig. 12 compares the second moments of the F,a(m) 
and F,7,(n) spectra. It is clearly shown that over 
1,500 cycles essentially no correlation exists between 
the u’ and v’ components, thus indicating that for n > 
1,500 cycles locally isotropic condition is approached. 


(7) Microscale Measurements 


Measurements of microscales of turbulence A, and 
A, (Fig. 13) standard 
placing two parallel hot-wires a distance Y or 


were accomplished by the 
method 
Z apart and measuring the correlation u’(0)u’(Y) or 
u’(O)u’(Z). From the curvature at Y = 0 or Z = Oof 
the correlation curve thus obtained, the microscales can 
be calculated. 

The distribution of \, was obtained by the differentia- 
tion method introduced by Townsend.'* In the lowest 
Reynolds Number case, because of the smal] differenti- 
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ated signal, no accurate measurements were possible. 
For R = 30,800, \, obtained from the spectrum meas- 
urements is also shown. From the fact that the corre- 
lation coefficient is a Fourier transform of the spectrum 
function, it follows that 


o?R 1 4 2 © 
( =) =—= oi n°F(n)dn 
OA sau Ae u- Jo 


Thus from the known distribution of the second moment 
of the spectrum function, \, can be easily calculated. 
This distribution is shown in Fig. 14 at different values 
of y/d. From Fig. 13 it is seen that the thus calculated 
\, values are in good agreement with those measured 
directly. 

Some preliminary measurements of an alternative 
method for obtaining A, are also indicated in Fig. 13. 
The method is described in detail in reference 14 and 
is based essentially en the fact that the number of zero 
values of a random function, such as wu’, is proportional 
to 1/A,. 

Fig. 11 indicates that the microscales decrease with 
increasing Reynolds Number as expected. Their distri- 
bution across the channel has a slight but definite 
maximum. 


CONCLUSION 


The present measurements indicate that a fully de- 
veloped channel flow can be realized experimentally. 
The semitheoretical predictions on the extent of the 
laminar sublayer and on the logarithmic distribution of 
the mean velocity are again proved to hold. 

The dimensionless velocity fluctuations, correlation 
coefficient k, and microscales decrease with increasing 
Reynolds Number. At the wall u’/u tends to approach 
a finite value independent of R. 

The measured microscales \,, \,, and A, show a con- 
sistent maximum at approximately y/d = 0.7. 

The shape of the u’ spectra obtained at different 
stations across the channel showed no appreciable dif- 
ference. 


Within certain limitation in accuracy described in 
the text, the high-frequency region of the «’ spectrum 
follows closely the n~—’ power law predicted by Heisen- 
berg. 

From the comparison of the u’ and shear spectra the 
existence of local isotropy in a turbulent channel flow is 
verified. 


REFERENCES 


1 Pai, Shih-I., Turbulent Flow Between Rotating Cylinders, 
N.A.C.A. T.N. No. 892, 1943. 

? Schubauer, G. B., and Klebanoff, P. S., Theory and Appli- 
cation of Hot Wire Instruments in the Investigation of Turbulent 
Boundary Layers, N.A.C.A. Wartime Report W-86. 

§’ Townsend, A. A., Measurements in the Turbulent Wake of a 
Cylinder, Proc. Roy. Soc., Section A, Vol. 190, p. 555, 1947. 

4 Corrsin, Stanley, Investigation of Flow in an Axially Sym- 
metrical Heated Jet of Air, N.A.C.A. Wartime Report W-94. 

5 Wattendorf, F. L., Investigation of Velocity Fluctuations in a 
Turbulent Flow, Journal of the Aeronautical Sciences, Vol. 3, 
No. 6, pp. 200-202, April, 1936. 

® Reichardt, H., Messungen Turbulenter Schwankungen, Na- 
turwissenschaften, Vol. 26, p. 404, 1938. 

7 Laufer, John, Investigation of Turbulent Flow in a Two- 
Dimensional Channel, Final Report, Contract NAw 5501, sub- 
mitted to the N.A.C.A. (in press). 

8 Goldstein, S., Modern Developments in Fluid Dynamics, 
Vol. I, p. 139; Oxford, 1943. 

® Fage, A., and Townend, H. C. H., An Examination of Turbu- 
lent Flow with an Ultramicroscope, Proc. Roy. Soc., Section A, 
Vol. 135, pp. 657-677, 1932. 

10 Heisenberg, W., Zur Statistischen Theorie der Turbulenz, 
Zeitschrift fur Physik, Vol. 124, 1948. 

11 Townsend, A. A., Local Isotropy in the Turbulent Wake of a 
Cylinder, Australian J. of Sci. Res. 1, No. 2, pp. 161-174, 1948. 

12 Corrsin, S., and Uberoi, M. S., Spectra and Diffusion in a 
Round Turbulent Jet, Final Report, Contract NAw 3505, sub- 
mitted to the N.A.C.A. 

13 Townsend, A. A., The Measurement of Double and Triple 
Correlation Derivatives in Isotropic Turbulence, Proc. Cambridge 
Phil. Soc., Vol. 34, p. 565, 1948. 

14 Liepmann, H. W., Laufer, J., and Liepmann, K., On the Spec- 
trum of Isotropic Turbulence, Final Report, Contract NAw 5632, 
submitted to the N.A.C.A. 











Crankshaft-Propeller Vibration Modes as 
Influenced by the Torsional Flexibility of 
the Engine Suspension 


B. M. FRAEYS pe VEUBEKE* 


Cwil Aeronautics, Brussels-Haren, Belgium 


ABSTRACT 


The present investigation is concerned with the fact that 
crankshaft-propeller torsional oscillations are coupled to the 
torsional oscillations of the engine as a whole. The basic phe- 
nomenon was discussed by Den Hartog and Butterfield! in the 
elementary case of a single radial engine with a single resonant 
frequency on the suspension. 

A general method of calculation is presented for long crank- 
shaft engines. It is an extension of Biot’s method,’ where, in 
addition to the usual dynamic moduli at both crankshaft ends, 
the solution contains the dynamic modulus of the engine sus- 
pension plus air frame with respect to the harmonic torque ex- 


erted by the crankcase. 


(1) INERTIA COUPLING BETWEEN CRANK ASSEMBLY 
AND CRANKCASE MOTION 


Bagi ¢, BE the angular coordinate of the xth crankpin 
with respect to the crankcase and y be the angular 
coordinate of the crankcase. The complete expression 
for the kinetic energy of the crank assembly is a quad- 
ratic, homogeneous expression in the angular velocities 
¢, and y 

T, = (1/2)b2I (bz) + ¥b2:M (bz) + (1/2)? V (dz) 


The coefficients are periodic functions, of which only 
the mean values J, /, and V will be considered. 

For instance, in the case of a single piston per crank- 
pin, these coefficients are : 


M = 1, + bR®{1 — (1/L)] 


1 R? 
.=M R? (1 x 
a +5 ( bs A 


H  J—bH(L — WH) 
p+ob I + — 


L? 
where 

R = radius of the crankpin 

I. = moment of inertia of crank section 
about crank axis 

p = mass of piston 

b, L, and JJ = respectively, mass, length, and dis- 
tance of center of gravity to crank 
axis for the connecting rod 

J = moment of inertia of connecting rod 


about the c.g. axis 


Received September 6, 1949. 
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It comes out that, contrary to V, the value of 1/ is the 
same for any location of the oscillation axis of the crank. 
case. This makes parts of J/, due to connecting rod 
and piston, additive in the case of several centrally 
connected pistons per crankpin (e.g., V-12 engines), 
The value of V is not given, since it depends on the 
oscillation axis and will be discussed later. 

The formulas are most easily obtained after a pre- 
liminary reduction to concentrated masses and cor 
rective inertias. Even the apparently involved case 
of the radial engine with angular eccentricities is 
tractable that way.' Referring to Tables 3 and 4 of 
this paper, where the formulas appropriate to the cal- 
culation of two coefficients Uy) and Vo are given, the fol- 
lowing correspondence exists: 


I = 2R°U ‘g, M = R°Vo Z 


(2) Eguations or Motion 


Fig. 1 shows an equivalent mechanical model of a 
long crankshaft engine with constant values of J and 
Let 6, denote the ampli- 
tude of vibration of the xth crankpin with respect to 
the crankcase and let a denote the amplitude of vibra- 
tion of the crankcase. 

The natural vibrations having the same frequency, 


M and a spring constant k. 


w, and the same phase angle, 6, + a, will be the ampli- 
tude of absolute motion of the xth crankpin. 

The equations of motion are then: for the «th crank- 
pin, 


—w7J6, — wMa = k(0,-1 + 6241 — 20,) (1) 
for the first crank, 


—w?lJ6; —_ w*Ma = k (02 ae 41) + My (2) 
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Fic. 1. Equivalent mechanical model of long crankshaft 


engine with constant inertia factors, including positive angle and 


torque conventions. 
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for the last crank, 


—wl0, — wMa = R(O,1 — 0,) — M, (3) 


for the crankcase, 


—w(I; + nV) nat w?V > 6, = Mo = M, te No _— 
1 
Nit N (4) 


where J, is the moment of inertia of the crankcase sep- 
arated from the propeller and blower gearing casings 
about the oscillation axis and where JN is the torque 
amplitude due to the engine suspension. 

(3) CoupLING DUE TO PROPELLER AND BLOWER 
GEARINGS 


The following linear and symmetrical relations hold 
between the torque amplitudes Ny and J/) and the 
corresponding absolute angular amplitudes a@ and 


6, + a: 


—M, => Ko(4, 4 a) 4. Goal 


— No /_ Go(8; + a) + Loa f (9) 


The coefficients are dynamic moduli whose expressions 
as functions of the frequency are to be set up in each 
particular case. With regard to the fixed frame prob- 
lem (a = 0), Ko is seen to correspond to Biot’s Ky.” 

The structures of the dynamic moduli in the case of 
a spur gear reductor are obtained as follows (Fig. 2): 
Let 


K, = —m,/ (0, + a) 


be the dynamic modulus at the end of the propeller 
shaft. The amplitudes 6, and 6,, relative to the gear 
casing and the torque amplitudes m, and my, due to 
teeth pressures, are related by the equations 


64 = TO»; My = ToMa 


(7) negative here) while 
No = (To = L)ma a wT a 


[ois the moment of inertia of the whole reductor group 
about the oscillation axis, the moment of inertia of the 
shaft groups about their own axis being ignored. 

On the other hand, 


Ma = Zaa(9a 4- = Za(9; + a) 
My = LZoalOa + “) = Zo0(9; + a) 


with Za = Zog 

The Z's are the dynamic moduli of the drive shaft 
“other end clamped.”’ Solving for the angular ampli- 
tudes, 


0a + a = (ma/Kaa) — 
A; 4. a= (Ma Kao) = 


(Mo Kw) 
(Mo/ Koo) 


The K’s are the dynamic moduli of the drive shaft 
“other end free’; they are simply related to the Z. 


ZwK aa = ZukKoa = ZaaK oo = ZaaZwo _ Zoa’ 
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Symbols used to set up dynamic moduli of spur gear 
type reductors. 





















































Fic. 4. Symbols used to set up dynamic moduli of epicyclic 
type reductor. 
These equations yield in succession 
64 + a = —(t) — L)ha — (107ma/K,) 
Ama = —Z 6; — [Za0 + = 1)Zaala 
where 


A= | ob (r0*Zaa K,) 
and, finally, 


Ko = Zw 1+ (to* Kea K,)] A- 
Go = ( . ] )\Zoa A 
Lo = ( ied l \*Zaa A 


1 _ wl, 
If the drive shaft has the simple structure of Fig. 3, 


Zaa = ky — wl 
Zu — Ro 
Zoo = Ro es wT 


If, furthermore, the inertias J, and /) are neglected, 


Laa = Za = Zw = Ro, Kae = 0 
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and there follows an exceedingly simple structure con- 
veniently referred to as the “‘purely elastic coupling”’: 


(1/hko) + (70°/ Ky) 
(7 = 1)Ko, ibe = (To = 1)*Ko —_ wl 


1/Ko 
Go 


One should be reminded that, owing to the flexi- 
bility of the propeller blades, A, and, hence, Ko will also 
be functions of the mean angular velocity of the crank- 
shaft.® 

The case of the epicyclic reduction gear is slightly 
more complicated. Fig. 4 shows elements similar to 
the spur gear reductor and similarly denoted. In 
addition to those, a crown wheel, of moment of inertia 
I,’, is connected to the gear casing through an elastic 
shaft of spring constant k, that allows in a crude fashion 
for the flexibility of the gear casing. Planet gears have 
a total moment of inertia about their own axis equal 
to T;. 

The relations between the different angles are 


4, = 


0, = 


— 19, + (I - Ho) 8. 
TOp ‘nal (To at 1)4, 


The same relations hold naturally when a is added to 
each relative angle. Using these relations in a kinetic 
energy equation, such as 


I,(6, + a)? = A(O, + a)? + BO, + a)? + C(Og + a)? 


values of A, B, and C are found that satisfy this equa- 
tion. It means that the total inertia J, of the planet 
gears may be transferred as follows: an amount 


A = (1 + wo) [1 + (uo/70) Zs 


goes to the crown wheel, whose fictitious total inertia 
becomes /,’ + A = /,; an amount 


5B = — pol (uo + t0)/(To = EY, 


goes to the planet carrier at the end of the propeller 
shaft; and an amount 


Cc = bo(1 + Ho) /To(T0 —_ 1) 


goes to the sun gear. 

In the expressions of the dynamic moduli that will 
follow it is understood that this transfer has taken place. 
There comes, in succession, 


i ies 


qn ee ( 1 
c= K. Ma K, Qa 


|(Ko/k — 1) cosp + 1 
\(K,/k — 1) cos mu + cos (n + 1)p 


The dynamic moduli Ko and K,, are functions of w and, 
hence, of u through Eq. (8). 
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64 ta= — |= : <i + Z| ma — (to — 1) ie 
with 
K, = ke — wl, 
and 
Amg = —ZwA — [Za + (to — 1)Zaalke/K,) a 
with 


A 


II 
oo 
N 
= 

pi, 
ia 
Fall 
| 
| 2 
eS is 
Ccmenomaid 


and, finally, 


tet eug§ Oe + SR a 
ie, aden ee 


Gy = (To = 1)Zoa(Re, i.) A- 
Lo = ( 1)°Z, (i) A-! — wil, — wl 
eo = (fo dd RK. W'Le K. Wt rp 


Similar formulas will hold for the other crankshaft 
end, eventually geared to blowers and (or) to exhaust 
turbines. 

M,, = KO, + Gra \ - 
M. +N, = Ch, > Leaf 6) 


(4) THE Frxep CRANKCASE PROBLEM 


Treated by Biot,” * this will correspond to the pres- 
ent problem where a@ is put equal to zero. The gen- 
eral equation of motion for the «th crankpin is satisfied 
by putting 

6.= Csin(xu +o) (x = 1,2...) (7) 
the frequency being related to the parameter yu by the 
equation 


ae , 
w = 2Vk/T sin (u/2) (8) 


Eqs. (2) and (3) are reduced to the following: 


[2 cos 4 — 1 a (Ko k)]6; ta 0: = 0 i] 9) 
[2 cos uw — 1 + (K,/R)]0, — On,-1 = OF 


and are used to fix the unknowns ¢ and uz. 

The substitution of the general solution [Eq. (7)] 
transforms them in a pair of linear and homogeneous 
equations in the unknowns sin ¢ and cos ¢. The 
condition of compatibility or “frequency equation” 
appears in the form of a determinant, 


(Ko/k — 1) sin w is 6 10 
(K,,/k — 1) sin nw + sin (nm + 1) yu! 


Instead of plotting the left-hand side as a trails 
cendental oscillating function of y to find its roots by 
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an interpolation process, Biot suggests constructing a 
complex quantity of which it is to be the imaginary 
part. It is useful, in view of the generalization of this 
process to the present problem, to draw attention to a 
convenient method for this construction. Any deter- 


minant 


' 4] = ad — bc 


is the imaginary part of either 


(a — 1b) (c + id) 


or 
—(a + ib) (c — id) 


As a result, the complex quantity is immediately re- 
solved in factors; in the case of Eq. (10), 


ele - 24 (Ko/k — er " x 


fe’? +. (K,/k — 1)e7*/?] = AoA ,e'™ + % Fn) (11) 


An equivalent statement to the conditior that the 
imaginary part should vanish consists in equating the 
argument of the complex quantity to a multiple of 7. 


Hence, 
nu + do + o, = multiple of 7 
where 
tan ¢ = (2k/Ko — 1) tan (u/2) ) } 
tan ¢, = (2k/K, — 1) tan (u/2)f =e) 


(5) LIMITATIONS OF THE SOLUTION 


It appears from Eq. (8) that the proposed solution 
yields only those natural frequencies that are inferior 


to the cutoff frequency 
/ 
o = 2VR/I 

of the mechanical filter represented by the long crank- 
shaft. This cut frequency is indeed the natural fre- 
quency of an elementary cell of that filter (Fig. 5). 

A corresponding process, yielding the natural fre- 
quencies, if any, above the cutoff frequency, would con- 
sist in putting 


6, = (—1)*C sinh (xu + 6) (x = 1,2...2) 
the general equation of motion being again satisfied, 
provided 


w = 2Vk/I cosh (u/2) 


When substituting the solution, note that 


| n 
~a? 2 0: = 2(cos wu — 1) 26 = 


The following additional notations are introduced: 


(cos nw — l)ny + sin @ — sin (u + ¢) — sin 





lp, 
, yw Vy 4 


td 4 





























6 40 


Fic. 5. Elementary cells of long crankshaft considered as a 
mechanical filter. The cut frequency is the natural frequency of 
a cell. 


However, it seems doubtful whether the consideration 
of these higher frequencies is worth while, especially 
since their accurate value will be impaired by the sim- 
plifications involved in the actual concentrated masses 
model. 

The same remarks concerning the limitation in fre- 
quency range will apply to the treatment of the general 


problem. 


(6) THE OSCILLATING FRAME PROBLEM 


Returning to the general equation [Eq. (1)], it is 
seen to admit of a solution 


6, = C[sin (xu + o) + vt 


4 2 ~ SI (13) 


C is again an arbitrary amplitude coefficient, and the 
frequency is still given by Eq. (8). The three param- 
eters u, ¢, and vy will follow from the equations govern- 
ing the motions or the end crankpins and of the crank- 
case. 

The investigation is limited to the case of 
elastic coupling’ on both ends of the crankshaft. The 
general case involves no special difficulties but is more 


“purely 


cumbersome. 
Substitution of Eqs. (13) and (5) in Eqs. (2) and (3) 
expressing the motions of the end crankpins yields 


(Ko, k— 1) sin (u + >) + sin @ — €)(Ko k)y = = uv] 


(K,/k — 1) sin (nu + o) + (14) 
sin[(m + 1)u +] — €n(K,/k)¥ = 0) 
where coupling coefficients appear 
ée = ™(I/M) — 1, @& =7,0/M) — 1 (15) 


Multiply Eq. (4) throughout by /(CkM)~". 


(mu + >) + 
sin [(7 + l)u + ¢] 


= I; + ny + To + ie 
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total moment of inertia of the engine about the oscillation axis; this makes a discussion about the exact value of y 
practically unnecessary, since the value actually included in an experimental measurement of F will, at least for g 
sufficient number of pistons, be independent of the angular position of the crankshaft and represent accurately 
the required mean 

K, = —N/a 
dynamic modulus representing the response of the engine mounting and air frame when submitted to a harmonic 
torque from the crankcase. This dynamic modulus is part of a more complicated expression that will be called the 


dynamic modulus of the suspension: 
h = —w*nl — (I/M)*[K, — w*F] (16) 
The coupling coefficients ¢ and ¢, are again put in evidence if Eqs. (14) multiplied, respectively, by 7o(// J) and 
T,(1/M) are substracted: 
é[sin (u + ¢) — sin d] + e,{sin (mu + o) — sin [(n + lu + o]} + (h/k)y = 0 (17) 


Eqs. (14) and (17) form a linear and homogeneous set in the unknowns sin ¢, cos ¢, and 7; hence, they yield a 


compatibility condition or frequency equation: 


\(Ko/k — 1) cosu + 1 (Ko/k — 1) sin uw —eKo/k | 
\(K,/k — 1) cos mu + cos (n + 1) (K,/k — 1) sinnu + sin (n+ 1)u —e,K,/k| = 0 
le(cos uw — 1) + €,[cos mu — cos (n + 1)u] eo sin + ¢,[sin mu — sin (x + 1)u] h/k | 
The development of this determinant is started through the elements of the last column, each determinantal 
cofactor being replaced by a complex quantity of which it is the imaginary part according to the process indicated 
for the fixed frame problem. Some terms, which are real, are conveniently dropped, and after some grouping of 


the others there comes: 


in h i(go + on) ° Ko - Pies F . . KoK,, 
er p Aodne i " + 2(1 — cos p) | &” + €,° — (l—e ™) (e? + e€,7) re — 


kb b 
Kok, ; 
Mu 
Eve, a @ real (18 
p? 
where 
Ko cos (u/2) K,, cos (u/2) 
Ap = and A, = 
k cos do k cos ¢, 


In contrast to the fixed frame problem, e’”” is no longer a factor of the whole expression; hence, the primitive idea 
of equating the total argument to a multiple of 7 will no longer give a satisfactory result, since the argument would 
depend upon circular functions of the angle mu, which means that it would go itself through many oscillations. 


However, it is possible to bring the vanishing con- i im +t im [X62 2 is 
Ae : ; ‘ ; b = cot an €)°Q, €,,7a¢) — aca,| (23 
dition of the imaginary part to a reduced form: 2 2 (eo"an F €n'o vat] 
3] 1 c= ( Tr . . . . 
sin (mu + ¥) + sin § 0 (19) The angle y is the argument of the bracket in Eq. (15), 


-. ; : which is written accordingly : 
where the angles y and ¢ are independent of my and vary ? 


smoothly. To this effect, introduce the dimensionless aKK 
flexibilities 0 n a@ Sil pb i(nu + y) — ip ’ 
i - : e — Nevene > real 
kk? 2 sin y 
2k 2k 2k 
ah , eos Qe So | (20) Finally, defining ¢ through one of expressions 
h a K,, ally, de £5 g e of expressions 
9 9) 9 
then put — Qere,rd Dern r Dene, r 
, sin { = —— siny = P — cosy = —== 
a ) V a2 bh? 
tan y = a/b (21) + (8 
where 


yields the reduced form [Eq. (19)] of the frequency 


a = a + a, — (&” + €,”)r (22) equation. 
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rg po Gn 
0 —90 —90 
1 —8.97 —68.05 
1.2 —4.915 —40.15 
1.3 —3.175 —7.74 
1.4 —1.58 +27 .38 
2.0 +6.14 +75.15 
5.0 +29.72 +86 .34 
7.0 +40. 23 +87 .47 
8.0 +44.69 +87.79 
9.0 +48. 52 +88 .06 
10 +51.815 +88 .26 
15 +64.13 +88 .86 
30 +82.785 +89. 44 
45 +101.838 +89 64 
60 +208 .055 +89. 74 
90 +256 .95 +89.85 
120 + 263.53 +89.91 
150 ; + 267.14 +89 .96 
180 +270 +90 


Fixed frame solution of the example shown in Fig. € 


TABLE 1 


Vi 


ny + po 


— 180 
—71 


—37 


—3. 


+34 
+93 
+ 146 


+169 .7 


+180 


+190.5 


+200 
+242 
+352 


+461. 
+657 .7 
+886. § 
+ 1,073 
+1,257.1 
+ 1,440 


4 


TABLE ¢ 
w,/27 = 88.63 Cycles per Sec. (u° = 8.05) 
“ vy” - nu+y +c 
0 — 180 0 — 180 
1.0 —77.26 —0.07 —71.3: 
1.5 +4908 —0.18 +57 .90 
2.0 +80. 44 —().09 +-92 35 
3.0 97.10 —(0.07 +115 
4.0 106.33 —().07 +130 
5.0 113.21 —0.08 +143 
7.0 134.36 —0.15 +176.% 
8.0 168.31 —0.73 +215 
9.0 308 . 44 —0.165 +362 
10.0 316.14 —0.065 +367 
15.0 332.46 —(0.01 +422 
30 351.94 +531 
15 371.44 +641 
TABLE 3 
w,/2r = 14.325 Cycles per Sec 
m y~ ig mu +y 
0 — 180 0 — 180 
l — 38.24 —5.40 —37 
1.2 -27 .745 — 16.59 +18.355 
1.3 +82 .438 — 22.39 +67.§ 
1.4 + 136.98 —18.72 +126 
1.5 +174.16 — 12.20 +170.96 
1.6 +197.16 —7.95 +198.: 
oe +212 .26 —5.49 +216.97 
1.8 + 222 895 —4.00 +229.7 
1.9 + 230.85 —3.05 +239 .2 
2.0 +237 .13 —2.41 +246 .72 
2.5 +256 .29 —1.02 +270 .27 


2.53 


26 
17 
40 
29 
04 
41 
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94 
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—Natural Frequencies 


“I 


31 


51.5 
66.! 


92.: 


121 
150 
Cut 


w/2r, 
Cycles 
per Sec. 


14.435 


337 .5 
548.6 
692.4 
912.2 
1,099 
1,220 


1,263 


frequency 


Natural Frequencies 

ng w/2r, 
Cycles 
per Sec. 


1.3) 14.435 
7.14 78.63 
8.96 98 63 
3] 337.5 


Natural Frequencies~ 


w/2r 


Cyc es 
p°? per Sec 
1.15 12.68 
1.44 15.87 
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Fic. 6. Inertias in lbs. in. sec.? and rigidities in lbs. in. per 


microrad. of typical V-12 engine with spur gear reductor. Fur- 
ther data used: 
I/M=1.3 

F = 700 lbs. in. sec.? 


There are two solutions to this equation 


nu + ~ + ¢ = even multiple of rl 
nu + w — ¢ = odd multiple of x f 


The process is thus reduced to the plotting of two 
smooth curves. To sum up: Plot the three expressions 
(20), the last two being already necessary in the fixed 
frame problem; then plot Eqs. (22) and (23), extract 
the angle y from Eq. (21), and, finally, extract ¢ from 
one of Eqs. (24). 


(25) 


(7) SprcrAL FREQUENCIES 


Some particular values of «4 are important to discuss. 
(a) IfK, = 0, then, 


tan y = cot (u/2)[—ao + AJ“! 


and sin ¢ = 0. 

The two curves have a common point if this condi- 
tion occurs for a particular value of uw (resonance at the 
rear dynamic modulus), or they are identical if the 
condition is permanent (rear crankshaft end free). 

The symmetrical case exists for resonant values of the 
front dynamic modulus. 

(b) If /# tends to zero, \ tends to infinity. 

Values of w for which this occurs are conveniently re- 
ferred to as resonant frequencies for the engine sus- 
pension, although the engine inertia and the inertia 
coupling factor J/.7 are involved in the expression for 


(8) SHAPE OF THI 


For each value of u ensuring the compatibility of Eqs. ( 


the shape of the natural mode in question. An obvious 


(14) and obtain linear and homogeneous equation in sit 
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Fic. 7. Two modal shapes of the engine in Fig. 6, when the 
resonant frequency of the suspension almost equals the second 
natural frequency of the same engine with fixed crankcase (88.63 
vs. 87.65 cycles per sec.). Frequency of first mode is 78.63 cycles 
per sec.; of second mode, 98.63 cycles per sec. Numbers refer 
to vibration amplitudes relative to the crankcase. To find abso- 
lute amplitudes, add values given for a. 


h. There comes 


tany = —(e? + €,”) cot (u/2) [eo?an, + €,?a0]~ 
sinf = — [2ee,/(€? + €,?)] sin 


(c) For antiresonant frequencies of the engine 


suspension, / tends to infinity, \ tends to zero, and 


(ao + a,) cot (u/2) 


tany = = tan (do + od» 


cot?(u/2) — aoa, 
snft = 0 
So we are naturally brought back to Biot’s fixed frame 


solution. 


=? VIBRATION MODES 

14) and (16), values of ¢ and y should be obtained to fix 
algebraic method is as follows: Eliminate y between Eqs. 
1¢and cos ¢. From this, tan ¢ is obtained, or, again, 4 


corresponding complex equation is set up—e.g., after some reductions 
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The argument x of the bracket is given by 


_ (€ + €,) tan (m/2)pu -|- (€ca@, + €,a0) tan (4/2) 
an x = ae eee ee 


(€) — €,) — (€9a@, — €,@) tan (u/2) tan (mp/2) 
and ¢ is determined by 
@+ [(n + 1)/2]u+x=0 
The determination chosen for ¢ is then used in one of Eqs. (14) to fix the corresponding value of y. From the 
first one, for instance, 
ey = cos (u/2) sin [¢ + (u/2)] — ao sin (u/2) cos [6 + (u/2)] 
(9) NUMERICAL EXAMPLE 


Fig. 5 shows the numerical values of inertias and duced. Flexibility of the propeller was not accounted 
rigidities pertaining to a V-12 engine with spur gear re- for, and the engine suspension was treated as a torque 
ductor. spring connecting the engine to a perfectly rigid air 

The interest being focused on the influence of the frame of extremely large inertia—hence, for practical 
engine suspension, some simplifications were intro- purposes, fixed in space. There is, accordingly, but 
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Fic. 8. Illustrating the determination of the two first natural frequencies of the engine in Fig. 6, when the resonant frequency 


of the suspension almost equals the first natural frequency of the same engine with fixed crankcase (14.325 vs. 14.445 cycles per 
sec.), 
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one resonant frequency for the engine suspension 
(h = 0), denoted by a,. 

Table 1 shows the results for a perfectly rigid sus- 
pension (fixed frame solution). For flexible suspension, 
a resonant frequency of 88.63 cycles per sec., very near 
to the second natural of the fixed frame case, was first 
investigated, since the frequency shift effect was ex- 
pected to be particularly large. 

Table 2 summarizes the results of calculations for 
this case. The angle y follows closely the sum @¢ + 
¢@, except in the neighborhood of the resonant fre- 
quency, where it gains 180° due to the added degree of 
freedom and then follows again the sum at 180° dis- 
tance (alternatively the Y curve might be notched as 
in reference 3). 

The angle ¢ is practically negligible. It shows two 
maxima: the first between the zeros of aj and a,; the 
second one near the resonant frequency of the suspen- 
sion. 

Figs. 6 and 7 show the calculated shapes of the two 
modes, one 10.3 per cent lower and one 12.5 per cent 
higher than the primitive second natural frequency. 
The other modes are practically undisturbed in fre- 
quency and shape. 


1950 


The case just discussed may be said to correspond to 
a ‘“‘hard”’ suspension. A case of “‘soft’’ suspension was 
next investigated by taking its resonant frequency at 
14.325 cycles per sec., which is near the first natural, 
As should be expected, only blower oscillations are 
affected with small frequency shift. 

The relatively large values of ¢, with the coalescence 
of the two previous maxima, are remarkable. 

Table 3 and Fig. 8 summarize the numerical results 

In all cases ¢ is exactly zero for u = 37.160, resonant 
frequency of the front dynamic modulus. 
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Solutions of the Lateral Equations of Motion 


(Continued from page 276) 


of-inertia terms gives rates of roll higher than those 
obtained by including the terms in the calculations. 
By defining a reduction factor k,, to account for side- 
slip and adverse yaw, such that 


(pb i gh = (Cy, Cy) “= Rr) 
the reduction in peak rolling velocity caused by sideslip 
and adverse yaw may be obtained as a function of lift 
coefficient. This has been done and is presented in 
Fig. 4, including and omitting the product-of-inertia 


terms. This also illustrates the effect of inclination of 


the longitudinal principal axis with respect to the 
flight path. 
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Plastic Collapse and Shakedown Theorems 
for Structures of Strain-Hardening Material’ 


B. G. NEALt 


Brown Un wersily 


ABSTRACT 


In recent years attempts have been made to develop methods 
of design for redundant structures of ductile material which are 
based upon the calculation of the load at which a structure ac- 
tually collapses as a result of excessive plastic deformation. Such 
methods of plastic or limit design may be used when it is known 
in advance that all the loads on the structure will be applied 
simultaneously. However, when a structure is subjected to 
several loads, each of which may vary independently of the 
other, a different type of failure must be considered. As the 
loads vary, a cycle of plastic deformations may be established, 
even though no possible combination of the loads could cause 
collapse. In this case, failure would occur either (1) by the frac- 
ture of some member in which plastic deformations occurred al- 
ternately in one sense and then in the other or (2) by the de- 
velopment of excessive plastic deformation as a result of a finite 
increase in the amount of plastic deformation occurring during 
each cycle of loading. It is therefore necessary to ensure that, 
as the loads vary, the structure attains a state of residual stress 
such that no further plastic flow can occur, no matter how often 
and in what order the loads are repeated. When this occurs, 
the structure is said to have shaken. down. 

Recent theoretical work on the problems of plastic collapse and 
shakedown has been restricted to the case of trusses and framed 
structures whose members exhibit the ideal-plastic type of stress- 
strain relation. In the present paper, some theorems concerning 
collapse and shakedown are established for trusses and framed 
structures whose members possess a strain-hardening type of 
characteristic. The most interesting cases, covered by the 
theorems, are those in which the strain hardening is limited, so 
that the stress-strain characteristic eventually becomes horizon- 
tal and plastic deformation can continue indefinitely under con- 
Stant stress. 

The method adopted for the proof of these theorems is to simu- 
late the strain-hardening characteristic by imagining a large 
number of bars of ideal-plastic material to be placed in parallel 
and constrained to have the same elongation. In this way 
theorems are established which are closely analogous to those that 
apply to structures of ideal-plastic material and which are in a 
suitable form to enable plastic collapse and shakedown loads to 


be calculated. 
INTRODUCTION 


M* METHODS OF STRUCTURAL DESIGN are based 
upon an elastic analysis of the structure, each 
member being so proportioned that the most unfavor- 
able combination of the external loads, when multiplied 
by a suitable safety factor, will just produce yield in 
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that member. However, it has long been known that 
a redundant structure is often by no means on the verge 
of collapse when yield occurs in one of its members.! 
If the external loads on a redundant structure are stead- 
ily increased, then, as one member yields, the structure 
may still be capable of bearing an increased load, those 
members that have not yielded then carrying a greater 
proportion of the total load. The purpose of what 
have been termed methods of plastic or limit design is 
to proportion the structure so that it would just col- 
lapse when subjected to the maximum specified loads 
multiplied by the chosen factor of safety. In this way 
the factor of safety takes on a definite meaning, whereas, 
if an elastic design method is adopted, the designer has 
no way of knowing whether application of the factored 
loads would bring the structure to a collapse state or 
would leave it capable of bearing a considerable increase 
in the loads. 

The earliest work on this subject appears to have 
been by Kist? and Griining.* Maier-Leibnitz‘ has 
carried out several important investigations on the be- 
havior of continuous beams loaded to collapse, and 
Baker has recently published a survey of researches 
carried out at Cambridge (England) which have re- 
sulted in the development of methods of plastic design 
for several types of frame structure.’ All these inves- 
tigations have involved the use of a principle that will 
be termed ‘the principle of plastic collapse.’’ This 
principle may be formulated as follows: ‘‘A set of ex- 
ternal loads will cause a structure to collapse if any in- 
crease of all the loads in proportion with one another 
would make it impossible to maintain equilibrium with- 
out violating the yield conditions.’ In the investiga- 
tions just described, the truth of this principle was usu- 
ally obvious, but for extremely complex structures 
this would no longer be true. However, a proof of 
this principle was recently given by Greenberg for 
the case of trusses and frames whose members exhibit 
the ideal-plastic relation between force and elongation 
or bending moment and curvature. This relation is 
of the form shown in Fig. 1 for the case of a truss mem- 
ber. 

It will be seen from this figure that the member is 
assumed to remain elastic unless the force has reached 
the tensile yield limit Y or the compressive yield limit 
Y’. When the force in the member is at either of these 
yield limits, the member can elongate or contract in- 
definitely under constant force. 
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Fic. 1. Force-elongation relation of ‘‘ideal plasticity’’ type. 

It often happens that a structure is subjected to 
several loads, each of which may vary between pre- 
scribed maximum and minimum limits irrespective of 
the values of the other loads at the same instant. In 
such cases, the structure should still be proportioned in 
such a way that no possible combination of the applied 
loads will cause collapse, but it is also necessary to en- 
sure that a cycle of plastic deformation is not estab- 
lished. Thus, the structure must be capable of finding 
its way to a state of residual stress such that further 
variations of load are supported in a purely elastic 
manner. When this occurs, the structure is said to 
shake down.’ A theorem governing the shakedown of 
trusses and frames whose members behave according 
to the ideal-plastic assumption was first conjectured 
by Bleich.’ This theorem may be stated as follows: 
“If any state of residual stress can be found for a struc- 
ture that enables all further variations of the external 
loads between their prescribed limits to be supported 
in a purely elastic manner, then the structure will shake 
down.” The actual state of residual stress in the struc- 
ture after it has shaken down will not necessarily be the 
state that has been found, and indeed it is possible that 
this latter state might be one that could never be 
reached by any possible loading program.’ This shake- 
down theorem was proved by Bleich’ for the case of 
trusses possessing not more than two redundant mem- 














ELONGATION 





Force-elongation relation of ‘“‘strain-hardening”’ type 
considered in paper. 


bers, and it was later established for the general case oj 
trusses with any number of redundant members py 
Melan.” A corresponding proof for the case of frame 
structures has been given by the present author.’ 

It should perhaps be remarked that, in many cases 
of structures subjected to repeated loads, other consid. 
erations may govern the design. For instance, in cases 
where the peak loads will be applied a large number of 
times, failure by fatigue must be guarded against, 
Design on the basis of shakedown would only be ra- 
tional if the peak loads were only expected to be ap. 
plied rather infrequently. 

The purpose of the present paper is to establish the 
principle of plastic collapse and the shakedown theorem 
for the more realistic case of structures whose members 
exhibit strain hardening. For this purpose the truss 
problem has been considered in detail, merely because 
of its greater mathematical simplicity. The corre- 
sponding results for the frame problem, which are 
stated at the end of the paper, may be derived in a pre- 
cisely similar manner. Before stating and _ proving 
the results, a detailed discussion of the assumed rela- 
tion between force and elongation for a truss member 
is given below. 


ForcE-ELONGATION CHARACTERISTIC 


Consideration will at first be given to trusses whose 
members are composed of material that exhibits a force- 
elongation characteristic of the type shown in Fig. 2. 
If a member is subjected to a steadily increasing tensile 
force S, the elongation A remains wholly elastic until 
the initial tensile yield limit Y is reached at the point 
a. In this region the force and elongation are con- 
nected by an elastic constant ¢ such that S = cA. 
Further increase in the force is accompanied by the de- 
velopment of plastic deformation, in addition to the 
elastic deformation. For instance, if the force is re- 
duced to zero after straining to the point ), a permanent 
deformation A, represented by Oc, will remain. The 
unloading process bc is purely elastic, so that the line bc 
is parallel to the original elastic line Oa. If the force ts 
further reduced from the point c so that the bar is com- 
pressed, plastic deformation will begin to develop at 
the point d. It will be assumed for the purpose of the 
present investigation that the elastic range bd remains 
the same, irrespective of the previous history of the 
deformation. Thus, if the loading program continues in 
such a way as to trace out the path defg, the elastic 
If the initial 
compressive yield limit for pure compression is }”, as 
shown, the change of force in any elastic range will 


range ef is assumed to be the same as bd. 


therefore be Y — Y’. 

If the force tends to a definite limiting value when 
the permanent deformation becomes large, the strain 
hardening will be said to be limited. This is the case 
illustrated in Fig. 2, where the extreme yield limit in 
tension is shown as Y, and in compression as Y,’. If, 
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STRUCTURES OF STRAIN 


1 the other hand, the force increases indefinitely as 


ol 
the permanent deformation increases, the strain hard- 
ening will be said to be unlimited, so that Y, = © and 
yi — OO 
yY, = . 


It is convenient to introduce a yield parameter }* 
such that at any stage of the loading program the ten- 
sile and compressive yield limits are Y + Y* and Y’ + 
y*, respectively. * is determined by the complete 
history of loading and is not a single-valued function 
of the permanent deformation. The elastic range of 
force, which is the difference between the tensile and 
compressive yield limits, is, of course, always Y — Y’ 
irrespective of the value of Y*. For the case of limited 
strain hardening, the requirement that for large pet- 
manent deformations the force tends to a limit imposes 
the following condition on }*: 
y,’-Y’'s Y*S Y,-Y (1) 
A further assumption is that, if yield is occurring, the 
increments of force and permanent deformation must 
always have the same sign. Thus, if S denotes the 
force and A denotes the permanent deformation, the 
increments S and A in these quantities must satisfy the 
condition 
SA 20 (2) 


the equality sign holding when the behavior is elastic 
or when the force is at one of the extreme yield limits. 

This latter assumption is, of course, in accordance 
with the actual behavior of engineering materials. 
However, the assumption that the elastic range of force 
remains constant is only true even approximately for 
certain materials. Investigations into the behavior of 
materials under the influence of cyclic loading were 
first carried out by Bauschinger,'! followed by Muir!” 
and Dalby.'* More recently, an extensive series of 
tests was reported by Templin and Sturm,'* and John- 
ston and Opila” have also carried out tests of this type. 
The results of such tests show that the effect of cyclic 
loading varies widely between different materials. In 
the case of structural mild steel, the sharply defined ini- 
tial yield point is removed after yielding has once taken 
place, and, if cyclic loading is continued, the force- 
elongation characteristic exhibits the well-known 
Bauschinger effect, thus reducing the elastic range by 
some 30-40 per cent. Templin and Sturm gave the re- 
sults of tests of tubular specimens of the aluminum alloy 
24S-T which were prestretched after heat-treatment 
and then tested in tension and compression. It was 
found that the elastic range actually increased by 
amounts up to 10 per cent as a result of the initial per- 
manent strain, which ranged from 0.2 to 3 per cent. 
On the other hand, tests by the same authors on solid 
circular specimens of 17S-T rod subjected to three com- 
plete loading cycles showed that the elastic range was 
reduced by about 10 per cent after the first yield and 
approximately constant. The 


thereafter remained 
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permanent strains produced in these tests varied up to 
0.6 per cent. Similar tests by Johnston and Opila 
on solid circular specimens of 27S-T rod subjected to a 
single loading cycle showed a corresponding reduction 
of about 25 per cent. 

Thus the theory of shakedown based on the assump- 
tion of a constant elastic range would require experi- 
mental confirmation before it could be used with confi- 
dence. Nevertheless, it is felt that the essential fea- 
tures of the behavior of structures subjected to repeated 
loading are described by the theory based on this sim- 
plifying assumption. 

It should be mentioned that, in quoting the above 
figures for the elastic range, hysteresis effects were 
ignored. This may be justified on the basis that a 
large number of applications of the various peak loads 
is not envisaged, since this would require the structure 
to be designed against failure by fatigue. Thus, the 
energy losses caused by hysteresis loops will not be ex- 
pected to be substantial in cases where shakedown has 


to be considered. 


THE PRINCIPLE OF PLASTIC COLLAPSE 


For proving the principle of plastic collapse, the 
following formulation is found to be most convenient. 
A set of external loads which would just cause a truss 
to collapse is characterized by the fact that, if every 
load were multiplied by any common factor greater 
than unity, it would be impossible to find any set of 
forces in the bars which would satisfy the conditions of 
both equilibrium and yield. 

Consider first a truss with » members, each of which 
is composed of ideal-plastic material, and let W be a 
typical external load. For convenience, imagine that 
all the loads maintain the same ratio with one another, 
so that W specifies the complete system of external 
loads. Let W, denote the set of external loads which 
causes collapse, and consider a load system ypW,, where 
uw is any positive number greater than unity. Then 
the principle of plastic collapse states that it is impos- 
sible to find any set of forces S,;* (¢ = 1, 2, ..., m) in the 
members of the truss which is in equilibrium with the 
external load system wl, and which also satisfies the 
yield conditions 

VY’ <S*< Y, 

In the case of a truss whose members exhibit strain 
hardening, the corresponding result that will be proved 
is stated in precisely the same way, except that the yield 
conditions become 

(Y,”. < S, S (¥,): (3) 


where (Y,’), and (Y,); are the extreme compressive and 
tensile yield limits for the 7th member, respectively. 
This result is clearly only of significance if the strain 
hardening is limited; in the case of unlimited strain 


hardening, Y,’ and Y, become —© and o, respec- 
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tively, so that the structure would never collapse in the 
sense that it would become incapable of bearing any 
increase of load. 

In establishing this result, the principle of virtual 
deformations will be used in the following form: Let 
S;/(t = 1, 2, ..., m) be any set of forces in the members 
which is in equilibrium under zero external load. 
Then, if A; (@ = 1, 2, ., n) is any set of compatible 
elongations in the members, 


n 
>, 5,4, = 0 (4) 
i=1 

Consider now a truss that has been loaded so that it 

has reached a state of collapse. Following Green- 
berg,® this is taken to imply that changes in the elonga- 
tions of the members take place while the external loads 
remain unaltered. Let S; be the changes of the forces 
in the members while in such a state of collapse. Then 
the corresponding changes in the elongations will be 
(S; /c;) + A,, where the S;/c; are the changes in the elas- 
tic elongations, and the A; are the changes in the per- 
manent elongations. Since the S; are in equilibrium 
under zero external load and (S; ci) + Aa being an 
actual set of elongations, must be compatible, substitu- 
tion in Eq. (4) gives 


Bs SAS; Ci) + Ail = 0 
t=1 


so that 


~ 


z. S22 q=- +o £4. (3 

i=1 i=l 
The left-hand side of Eq. (5) is clearly positive defi- 
nite, whereas Eq. (2) shows that the right-hand side 
It follows that both sides of the 
Thus, 


is negative definite. 
equation must be zero, so that all the S; are zero. 
when the truss is in a state of collapse, the force in 
every member remains constant as the deformation 
under constant external load continues, and the changes 
of deformation which occur are entirely plastic. For 
those members that continue to deform in the collapse 
state, the force must be equal to one or other of the 
extreme yield limits, Y,, or Y,’, since otherwise the force 
in these members would also change in accordance 
with the strain-hardening property of the material. 

Suppose now that, contrary to the principle of plastic 
collapse, a set of forces S;* could be found which is in 
equilibrium with the load system ulV, and does not 
violate the extreme yield limits for any member, so 
that 


ey? < oe < 2 AP 


If any such set of forces actually existed, there would 
also exist a set of forces S;*/u which would be in 
equilibrium with the loads IV,, and it is clear that this 
set of forces would not violate the yield conditions. 
Thus, the set of forces [S; — (S;*/u)| would be in 
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equilibrium under zero external load, S; being the actual 
set of forces in the structure at collapse. Using this 
set of forces in Eq. (4), together with the compatible 
set of deformations A,, it follows that 


n 
dS: — (S;*/p) JA; = 0 (6 
c= 

Now, if for any member A; > 0, the actual force S; in 
the member must be the extreme tensile yield limit 

(Y,)i, since otherwise this force would change as the 

plastic deformation changed. Thus, if A; > 0, S, = 

(Y,); and S;*/n < (Y,);, and, similarly, if Ai < 0, 

S; = (Y,’); and S,*/pn > (Y,’):. 


It follows that, for each member, 
[S; “ie (S;* wu) Ai 2 0 


Comparing this result with Eq. (6), it is seen that, for 
each member, 


[S; — (S:*/u) JA; = 0 (7 


Now, when the truss is in a state of collapse, there 
must be at least one member, say the kth member, in 
which permanent deformation is occurring. For this 
member, A; ¥ 0, so that S, = S,*/u and either 


Si _ (Yaa S,* < (Yp)x 
or 


Si _ (Vp )is S.* 2 (Fy be 


Pg 


In either case it is at once seen that uw < 1, thus vio- 
lating the original hypothesis that u > 1. This im- 
plies that no set of forces S,;* exists, so that a set of ex- 
ternal loads which causes collapse is characterized by 
the fact that these loads cannot be increased in propor- 
tion with one another while still permitting the condi- 
tions of equilibrium and yield to be satisfied simultane- 


ously. 


SHAKEDOWN THEOREM 


Consider a truss composed of » members which is 
subjected to several external loads, each of which can 
between prescribed maximum and minimum 
values. At any stage in the loading program, let 5; be 
the force in the ith member, and let B; denote the force 


that would be imposed on this member if the entire truss 


vary 


were behaving elastically under the same external loads. 
It will be convenient to define the residual force in the 
ith member as 


R, = S; — B; (d 


The residual force defined in this way is the force that 
would arise in the ith bar upon removal of the external 
loads, provided that the truss behaved elastically dur- 
ing the unloading process. It is, of course, possible 
that further yielding might take place while the loads 
were being removed, in which case the R; would not 
have this simple physical significance. 
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Assuming the truss to behave elastically throughout, 
the force in any given member corresponding to unit 
yalue of any of the external loads can be calculated 
readily by standard methods. By suitably combining 
these unit elastic solutions with either the maximum 
or minimum values of the external loads, the greatest 
and least possible elastic forces in any member due to 
variations of all the external loads between their pre- 
scribed limits can then be calculated, using the princi- 
ple of superposition. For example, if, because of unit 
positive value of the external load W,, a positive (ten- 
sjle) force tf, is caused in the /th member of the truss, 
the contribution to the maximum elastic force in this 
member would be ¢,W,.""", where W,” is the pre- 
scribed maximum value of the external load VW’. 

Let B;”"™: and B,”*" denote the maximum and mini- 
mum values of the elastic force in the 7th member of the 
truss as calculated inthis way. Then the shakedown 
theorem for trusses whose members are composed of 
ideal-plastic material may be formulated in the follow- 
ing way: If any particular set of residual forces R; 
can be found which satisfies the inequalities 


mar. RD. < ’ 
B; = is R; ~ Y; \ (9) 


Bi". 4. R, > y 


for all members of the truss, shakedown will occur. 
Y, and Y,’ are, of course, the tensile and compressive 
yield limits, respectively, for the 7th bar of the truss, 
and the R,; are required to satisfy the conditions of 
equilibrium under no external load. It should perhaps 
be pointed out that the actual set of residual forces in 
the truss after it has shaken down may very well differ 
from the R;, and indeed it is possible that a particular 
set of residual forces R; which has been found to satisfy 
the inequalities (9) could never be reached by any load- 
ing program whatsoever.’ 

Theshakedown theorem in this form cannot be applied 
to the case of trusses composed of strain-hardening ma- 
terial, for, as the various members yield, their yield 
limits change, as previously discussed. However, the 
following theorem can be established for trusses of this 
type: If any set of residual forces R, can be found 
which satisfies the inequalities 


or = a R, S$ Y; 5 i Y,* | (10) 
By" +R, 2 Y, + ¥,°9 
where 
(Y,’), -— Y,’ < ¥* < (Y,): — Y, (11) 


shakedown will occur. Here the Y;* are any set of 
yield parameters, each of which satisfies the inequalities 
(11). 

If the members of the truss exhibit unlimited strain 
hardening, so that Y, and Y,’ become © and —o, 
respectively, there ceases to be any restriction on Y,* 
through the inequalities (11). The inequalities (10) 
May then be written in the form 
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Fic. 3. Linear strain-hardening characteristic assumed by 
Melan. 


Bm. + R, — ¥, < Y,* < By + R, _ ry 


and, remembering that Y,* is unrestricted, the condi- 
tion for shakedown reduces to 


, >? 4 
B,*- — B,**. < Y; _ y, (12) 


This result was established by Melan" for the special 
case of linear strain hardening, where the relation be- 
tween force and elongation for any member of the truss 
is as Shown in Fig. 3. 

To prove the theorem that shakedown will occur if 
the inequalities (10) and (11) can be satisfied, it will be 
convenient to imagine the truss to be replaced by a 
model in which each member of the truss is replaced 
by an assemblage of parallel bars of ideal-plastic ma- 
terial. It will be shown that such assemblages can be 
made to reproduce the essential features of the force- 
elongation characteristics of members composed of 
strain-hardening material when subjected to any pro- 
gram of cyclic loading. Thus, if it can be shown that 
the model will shake down, it follows that the actual 
truss would also shake down under the same loading 
conditions. It should be mentioned that this device 
has been used by Bohnenblust and Duwez" for a dif- 


ferent purpose. 


A MODEL FOR STRAIN HARDENING 


Consider an assemblage of a large number of parallel 
bars, constrained at their ends to have the same elonga- 
tion. Each of these bars is composed of ideal-plastic 
material, and they each have the same elastic constant 
but different yield limits. It will be shown that the 
overall force-elongation characteristic of such an as- 
semblage can be made to reproduce the essential fea- 
tures of the force-elongation characteristic of a strain- 
hardening member under cyclic loading. For the 
sake of simplicity, a detailed discussion will first be 
given for the case of two parallel bars. 

Consider two parallel bars, each possessing the same 
elastic constant, and suppose that bar | yields under a 
smaller force in both tension and compression. The 
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yield limits in tension are denoted by ; and ye, and in 
compression y;’ and y2’, so that y; < ye and y,’ > yo’. 
The resultant force is, of course, the sum of the indi- 
vidual forces in the two bars, and the elongation is the 
same for each bar. Fig. 4 shows the behavior of each 
of the bars as the elongation varies and also the result- 
ant force-elongation behavior. The resultant force is 
plotted to half-scale in this figure, so that for a given 
elongation the resultant force is obtained by taking 
the mean between the values plotted for the individual 
forces. The variation of the resultant force for a par- 
ticular variation of the elongation is shown by the 
dotted line Oabcdefg. Along Oa, both bars remain elas- 
tic, bar 1 reaching its tensile yield limit y, ata. As the 
elongation continues to increase, bar 2 remains elastic 
until the point 5 is reached, and then along bc both the 
bars are yielding. At c, the elongation is reduced 
until the point f is reached, when it is then increased 
again. A detailed description of the behavior is not 
necessary, but a few points of particular interest are 
brought out. 

(1) Irrespective of the previous loading history, bar 
2 can never yield unless bar 1 is also yielding. 

(2) The elastic range, cd, fg, etc., remains constant 
at the value 2(y, — 9’). 

(3) The length of the ‘“‘elbow,”’ 
sum of the lengths of the two elbows, ab and 1, oc- 
curring at first yield in tension and compression, re- 


de, is equal to the 


spectively. 
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Behavior of system of two parallel bars having different 
yield limits. 


Fic. 4. 


These results may be generalized to cover the case 
of m parallel bars. Thus, if bar 1 yields under the 
smallest force in both tension and compression, yield 
can only occur in any of the bars if bar 1 is yielding. 
This result is of particular significance in the work that 
follows. 
value irrespective of the previous loading history. By 


The elastic range also preserves the same 


letting m tend to infinity and by arranging a suitable 
distribution of yield limits among the bars, the behavior 
of an actual strain-hardening material could be approxi- 
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mated in this way, although it should be noted that the 
shape of elbows of the type corresponding to de in Fig, 
4 would be determined completely by the shape of the 
elbows corresponding to ab and fi in Fig. 4. This 
would not, of course, be the case for real materials, byt 
this discrepancy does not invalidate the results that 
will be presented later. 

In order to obtain some results that will be required 
for the proof of the shakedown theorem, consider the 
case of m parallel bars in more detail. Suppose that 
the assemblage corresponds to the 7th member in a truss, 
which possesses an elastic constant c;. Since each bar 
of the assemblage possesses the same elastic constant, 
the value of this constant must be c;/m. If a tensile 
force is applied to the assemblage of bars, starting from 
a stress-free state, yield will first occur in bar 1. When 
this takes place, the elongation of bar 1, and thus of 
every bar, will be my,/c; Thus, the force required to 
produce initial tensile yield is my, so that 

my, = Y, \ 
my,’ = Y;,’ f (18) 

If the forces in the individual bars of the assemblage 
are denoted by s,(p = 1, 2, ., m), equilibrium de- 
mands that 


m 
S,=> 5, (14) 
p=1 


Since each of the bars of the assemblage has the same 
elastic constant, it is clear that 


m 


B, = > b, = mb, (15 
p=1 


where b, is the elastic force in the pth bar of the assem- 
blage. 

From Eqs. (8), (14), and (15), it follows at once that 

m 
R; = 8 lp (16 
p=1 
where 7, is the residual force in the pth bar of the assem- 
blage, defined as s, — dy. 

The residual force R; is the force in the ith member of 
the truss when the external loads on the truss are re- 
moved, assuming that the truss behaves elastically dur- 
ing the unloading process. If this member were then 
removed from the truss, the force on the member would 
be reduced to zero, but the equivalent assemblage 
would possess a set of residual forces which will be de- 
noted by 7,’.. Remembering that unloading processes 
are assumed to be elastic, it will be seen that, if the as- 
semblage were replaced in the unloaded truss, the force 
R, would be shared equally between the m bars, so that 
the force in the pth bar would ber,’ + (R;/m). Let this 
force, which is the residual force in the pth bar of the 
assemblage when the truss is free from load, be denoted 


by r,. Then, 


tp = Tp + (Ri/m) (17) 
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so that 


m m m 
> % = > [r>’ + (Ri/m)) = > 7, + BR 
p=1 p=l p=1 


Since the 7,’ are in equilibrium with zero load on the as- 
m 
semblage, >. 7,’ = 0, so that this result agrees with 
. p=1 
Eq. (16). 
Suppose now that, starting from the state when free 


from external load, with residual forces R;, the truss is 
loaded so that the 7th member sustains an additional 
elastic force P;, which just causes yield in tension. 
Equating the total force on the 7th member to the ten- 
sile yield limit, it is found that 


Pi + & = ¥,+ ¥,° (18) 


The additional elastic force P; is shared equally among 
the bars of the assenrblage. Since yield in tension is 
just occurring, bar 1 must just be at its yield limit. 
Equating the total force on bar 1 to its tensile yield 
limit, it is found that 


(P,/m) +n =n (19) 


Eliminating P; between Eqs. (18) and (19) and using 
Eqs. (13), it then follows that 


mr, = R, — Y,* (20) 


Considering the pth bar of the assemblage when yield 
in tension is just occurring, it is noted that the force 
in this bar must be less than, or equal to, its tensile 
yield limit, so that 


(P,/m) + tr < (21) 


Combining this result with Eq. (18), it then follows 
that 


mr, < my, + Ri — Yi — Yi* (22) 
By a similar argument it can readily be shown that 
, ry r o« 
mr, 2 my,’ + R; — Y;' — Yi* (23) 


The significance of the results contained in Eq. (20) 
and the inequalities (22) and (23) is that to any pair of 
values of R; and Y;,* there always corresponds at least 
one set of values of r, which satisfies these relations, 
provided, of course, that the Y;,* satisfy the inequalities 
(1); for it is clear that it is possible to attain any such 
value of Y,* by at least one loading program, and thus 
with any value of Y,* there can be associated a set of 
residual forces r,’ which would appear in the assemblage 
when free from load. When the assemblage is imag- 
ined to replace the ith member of the truss in which the 
residual force is R,, the force on the assemblage becomes 
R,, and thus each member of the assemblage is sub- 
jected to a further force R;/m. Thus, for any pair 
of values of Y,;* and R, there exists at least one 
set of residual forces 7, in the assemblage such that 
t = 1,’ + (R,/m) satisfies the relations (20), (22), ard 
(23). 
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PROOF OF THE SHAKEDOWN THEOREM 


In order to establish the shakedown theorem for 
trusses whose members exhibit strain hardening, each 
member of the truss is imagined to be replaced by an 
equivalent assemblage of ideal-plastic bars, as just dis- 
cussed. If it can be shown that the model thus ob- 
tained would shake down under the same set of ex- 
ternal loads, it follows that the actual truss would also 
shake down, for each member of the truss possesses a 
force-elongation characteristic that is identical with 
that of the corresponding assemblage. The advantage 
of this procedure is that the model is a truss in which 
every member is composed of ideal-plastic material, so 
that the shakedown theorem for such trusses as em- 
bodied in the inequalities (9) is applicable. 

Suppose that for a particular truss it has been found 
possible to satisfy the inequalities (10) and (11). From 
Eq. (20) it is known that to any pair of values of Y;,* 
and R; there corresponds a definite value of the residual 
force 7; in bar 1 of the assemblage. Thus, from the 
inequalities (10), 


B,"™ + mn < Y; 1 _ 
By. + mr, > Vif (24) 
where fr; is the residual force in bar 1 of the assemblage 
which corresponds to the chosen values of R, and Y,*. 
Dividing these inequalities throughout by m and using 
Eqs. (13), it is then seen that 


(B,"™-/m) +n Sy i] 


. min. / - oa? (25) 
(Byn™"-/m) + 7% 2 ni’ 


Now (B,”"-/m) is the maximum elastic force 6,” 
appearing in any bar of the assemblage, since each bar 
possesses the same elastic constant, and (B,""/m) 
admits of a similar interpretation. Thus, 

bh" +n in | 26) 
, { «V0 
5, **- a ial > yi’ § y 

It is also possible to combine the inequalities (22) and 

(23) with (10), giving 
Bywv™ £ Y¥,+ Yi* — Ri <& my, — mr, 


By" 2 Y, + Y,* — R; 2 my,’ — mr, 


so that 


_-" + lp > Vp \ (27) 


~ + ry v vy’ S \< 


where 7, is a residual force in bar p of the assemblage 
which corresponds to the chosen values of R, and Y,*. 

Thus, if each member of the truss is considered to be 
replaced by an assemblage of bars which behaves in an 
equivalent manner under cyclic loading, the foregoing 
argument shows that for the model comprising the as- 
semblages the inequalities (26) and (27) can be satis- 
fied, provided that the inequalities (10) and (11) can 
be satisfied. However, as already pointed out, the 
shakedown theorem for ideal-plastic trusses is applicable 
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to the model, and this theorem states that, if the in- 
equalities (26) and (27) can be satisfied, the model will 
shake down. Thus, it follows that, if the inequalities 
(10) and (11) can be satisfied, the actual truss will shake 


down. 


FRAME STRUCTURES 


Consider now the case of frame structures that carry 
load by virtue of the resistance to flexure of the mem- 
It will be assumed that the bending moment- 
at any cross section is similar 


bers. 
curvature characteristic 
in form to the force-elongation characteristic shown for 
a truss member in Fig. 2. Thus, for positive curvature, 
there will exist an initial yield moment My at which 
plastic curvature begins to develop and a fully plastic 
M, developed when the curvature becomes 
the corre- 


moment 
For negative curvature, 
M,’, respectively. 


extremely large. 
sponding quantities are My’ 
The elastic range of bending moment is assumed to re- 
My’, and a yield moment parani- 


and 


main at the value My — 
eter 1/* may be defined such that after any loading 
program the yield moment for positive curvature is 
My + M* and for negative curvature is My’ + M*. 
This parameter must, of course, satisfy the conditions 


M,’ — My’ < M* < M, — My (28) 


The analogy that exists between trusses and frame 
structures that carry loads by virtue of the resistance 
to flexure of the members has been discussed by the 
author in a previous paper.’ If it is assumed that the 
frame is subjected only to concentrated loads, the bend- 
ing moment diagram for any member will consist of a 
series of straight-line segments that terminate at what 
are termed the joints in the structure. Since a linear 
variation of bending moment along a beam implies 
that the shear force is constant, a joint is seen to occur 
at any section where a transverse force is applied to a 
member. Thus, joints occur at any section where an 
external load or couple is applied to a member or where 
the member rests on a support or is connected to an- 
other member of the frame. 

It is then assumed that yield cannot occur in the span 
of a member between two joints unless yield is occur- 
ring at at least one of the joints. In many structures 
the cross section of the members between the joints 
remains constant, so that this assumption is obviously 
justifiable. If the member is of variable cross section, 
it may still be possible to retain this assumption.’ 

With these assumptions, the state of stress in the en- 


tire frame is specified by the values of the bending 
moments at the joints, and, if it can be shown that yield 
is not occurring at the joints, it is known that yield is 
Thus, if there 
are m joints in a frame at which the bending moment 


joint 


not occurring anywhere in the frame. 


is nonzero—thus excluding, for example, the 
which occurs at the end of a beam resting on a simple 


support—it would be expected that the principle of 


‘ 
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plastic collapse could be established in the following 
form, exactly analogous to the case of a truss: A set 
of external loads which would just cause a frame to 
collapse is characterized by the fact that, if every loaq 
were multiplied by any common factor greater than 
unity, it would be impossible to find any set of bending 
moments at the joints which would be in equilibrium 
with the increased loads without violating the yield 


conditions 
(M,’); < M, g (M,): (29) 


This result can readily be established in precisely the 
same manner employed for the case of a truss, and for 
the sake of brevity the proof is omitted. 

It would also be expected that a shakedown theorem 
could be established in the following form: Let my" 
and ,;""" be the maximum and minimum bending 
moments occurring at the ith joint of the frame if 
wholly elastic behavior is assumed. Then, if any par- 
ticular set of residual moments 7; can be found which 
satisfies the inequalities 


My"? + m, < (My); + M,* \ ' 
min = , 7 x (30) 
Mi, +m, 2 (My’), + M, f 


where 


(M,'); — (My’); & M* < (M,);i — (My); (31) 


the structure will shake down. This result can also be 
established, as in the case of a truss, with the aid of a 
model. Details of the proof will be omitted, but it is of 
interest to note the kind of model that can be used for 
this purpose. 

Consider first the behavior of a pair of parallel rods 
of strain-hardening material whose axes are separated 
by a distance 2h, as shown in end view in Fig. 5. The 
rods are assumed to be of small diameter compared with 
the distance 4, and they are imagined to be connected 
by a thin web, indicated by the dotted line in Fig. 5, 
which transmits the shear force necessary to cause a 
linear variation of bending moment along the beam. 
The behavior of tension is assumed to be 
identical with the compres- 


rod 1 in 
behavior of rod 1’ in 


sion, and vice versa. Thus, a bending moment J/ ap- 
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plied about the axis AA is resisted equally by either 
rod, and, if the force in rod | is denoted by S, 


M = 2hs (32) 


assuming that the web does not contribute to the bend- 
ing moment. Furthermore, it is readily seen that, if 
the curvature is x, the longitudinal strain in either rod 
is given by 

c=e/h (33) 
Eqs. (32) and (33) show that the bending moment-cur- 
vature characteristic is precisely similar in form to the 
force-elongation characteristic for either rod. 

To obtain the actual model, let each of the strain- 
hardening rods be replaced by an equivalent assem- 
blage of parallel rods of ideal-plastic material, each 
possessing the same elastic constant but having different 
yield limits. Such an arrangement is illustrated in 
Fig. 6, in which, for the sake of simplicity, the case of 
only two pairs of ideal-plastic rods is illustrated. In 
order to preserve symmetry, the ideal-plastic rods are 
arranged in pairs about the principal axis BB, so that 
when yield occurs there is no resultant bending moment 
about this axis. Since each rod must have the same 
longitudinal strain at every cross section, a large num- 
ber of rigid transverse rods are arranged as shown in 
the plan view in Fig. 6. 

It is clear that in this way a model is obtained whose 
bending moment-curvature characteristic at any cross 
section is precisely similar in form to the force-elonga- 
tion characteristic of the assemblage of bars in either 
flange. Thus, results equivalent to those obtained for 
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the truss model can be written down at once, enabling 
the shakedown theorem contained in the inequalities 
(30) and (31) to be established. 


CONCLUSIONS 


The shakedown and collapse theorems established 
in this paper for the case of a truss whose members 
exhibit strain hardening are a little more realistic than 
the corresponding theorems for ideal-plastic trusses. 
The only material that approximates to the ideal-plastic 
assumption is mild steel, and a compression member of 
ideal-plastic material would buckle as soon as the com- 
pressive yield limit was reached and thereafter con- 
tinue to deform axially at a greatly reduced load. This 
was pointed out by Symonds and Prager,’ who stated 
that their object in considering the ideal-plastic truss 
was to present basic ideas in the simplest possible way, 
rather than to establish results of immediate practical 
value. A similar philosophy forms the justification 
for that part of the present paper which deals with 
strain-hardening trusses, for a compression member of 
strain-hardening material would buckle before the ex- 
treme compressive yield limit was reached, and axial 
deformation would then continue at a lower load.'* 

A further point should be mentioned in connection 
with the behavior of compression members. Consider, 
for the sake of simplicity, the ideal case of a slender 
pin-ended strut that is perfectly straight and loaded 
axially. Below the Euler critical load, the axial de- 
formation is proportional to the load. When the Euler 
load is reached, the lateral deflection and, therefore, 
the axial deformation become large, while the load 
remains constant. Thus, the force-elongation relation 
for steadily increasing deformation is precisely of ideal- 
plastic form. However, if the strut is sufficiently slen- 
der, this behavior would be purely elastic, so that, if the 
deformation were reduced, the load would remain con- 
stant until the central deflection became zero, in con- 
trast with the ideal-plastic assumption of elastic un- 
loading. Such an effect would also be noticeable in 
struts that had partially yielded during buckling. 

These effects restrict the validity of the theorems 
proved for trusses to those cases in which the compres- 
sion members are braced against lateral deflection. 
The corresponding results for frames are not restricted 
in such a direct manner, in the sense that there is no 
distinction between positive and negative curvature. 
However, in order for the shakedown and collapse 
theorems to be applicable, it would still be necessary to 
ensure that no member in a frame would buckle under 
axial load and, also, that no member would be suf- 
ficiently thin and deep to buckle laterally when bent 
about its stronger principal axis. With these restric- 
tions, the principle of plastic collapse for frames is 
suitable for direct application, since beams of structural 
mild steel and, also, the light alloys exhibit bending 
moment-curvature characteristics of the type shown in 
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Fig. 2, except that the elastic range is not usually con- 
stant. However, this assumption did not enter into 
the proof of the principle of plastic collapse. The 
proof of the shakedown theorem, on the other hand, 
rested directly on this assumption, so that this result 
could not be used without the aid of experimental work. 

Unfortunately, it was not possible to include any 
numerical examples to illustrate the application of these 
theorems because of limitations of space. However, it 
will be seen that in each case the investigation would be 
concerned with finding whether any particular solution 
could be found for a set of linear inequalities. A gen- 
eral method for investigating problems of this type has 
been given by Dines,'® and the application of his 
method to particular problems is being examined. 

While the theorems presented in this paper enable 
the load-carrying capacity of structures to be examined, 
an important question has been left unanswered. It 
might happen that a structure, while still theoretically 
capable of carrying further load, would become prac- 
tically useless because of excessive deformations of some 
of its members. Intuitively, it seems clear that, while 
yielded members are still contained by a framework of 
elastic members which in itself is capable of carrying 
load, excessive deformations could not occur. It 
would, however, be desirable to be able to calculate 
an upper limit on the deformations of any given mem- 
ber under these circumstances. 
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Leading-Edge Singularities in 
Thin-Airfoil Theory 


ROBERT T. JONES* 
Ames Aeronautical Laboratory, N.A.C.A. 


SUMMARY 


If the thin-airfoil theory is applied to an airfoil having a rounded 
leading edge, a certain error will arise in the determination of the 
pressure distribution around the nose. It is shown that the 
evaluation of the drag of such a blunt-nosed airfoil by the thin- 
airfoil theory requires the addition of a ‘‘leading-edge force,”’ 
analogous to the well-known leading-edge thrust of the lifting air 
foil. The method of calculation is illustrated by application to 
1) the Joukowski airfoil in subsonic flow and (2) the thin elliptic 
cone in supersonic flow. “The paper concludes with a general 
formula for the edge force which is applicable to a variety of 


wing forms. 
INTRODUCTION 
taken to 


evaluate the effects of singularities that appear in the 
One example of such an effect is the 


- THE APPLICATION of thin-airfoil theory to air- 
foils of be 


finite thickness, care must 


thin-airfoil flows. 
finite force arising from the singularity at the leading 
edge of a lifting airfoil—that is, the well-known “‘lead- 
ing-edge thrust.” 

Another example is the finite force arising from the 
singularity in the flow field of a nonlifting airfoil hav- 
ing a blunt leading edge. The latter effect does not 
appear to have been discussed previously and is of im- 
portance at supersonic speeds where it is desired to 
calculate the arising from the thickness 
of the airfoil. 


“wave drag 


THIN JOUKOWSKI AIRFOIL IN SUBSONIC FLOW 


A simple example of the leading-edge force is pro- 
vided by the two-dimensional flow over an airfoil of the 
Fig. | shows such an airfoil and the 
According to the thin-airfoil 


Joukowski type. 
coordinate axes used. 
theory, the conjugate u-iw of the perturbation velocity 
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Fic. 1. Thin Joukowski airfoil. 
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vector will be represented by an analytic function of the 


complex variable § = x + iz; that is, 


u — iw = f(é) 
The function f is chosen to satisfy the desired boundary 
conditions around the slit —1 < « < +1 and to make 
u — tw vanish at infinity. 
By a process of trial, it was found that for the example 


chosen 
u— w= 
, l — é J 
Vt +ave ~1-8) +1 (1) 
Ve — 1 
Near the chord line £ = x + ot and 
u/Vt > 1— 2x 
WwW l aad X 
> + ( — 2V1 — x? (2) 
Vt V1l-x 


The ordinate of the airfoil surface, obtained by integrat- 
ing w/V = dz/dx is 
(3) 


z= 1 —x)V1 — x? 


The constant ¢ is thus the ordinate of the airfoil at the 


50 per cent chord station (x = 0). The pressure co- 
efficient is given by 

Ap/(1/2)pV? = —2(u/V) (4) 
and, hence, is proportional to 1 — 2x in the interval 


—1<x< +1. At the ends of the interval the terms 
under the radical change sign and the expression for 
u — tw become a pure real number, so that w = 0 
along the streamline ahead of and behind the airfoil. 
Fig. 2 shows the pressure distribution along this stream- 
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Fic. 2. Pressure distribution for Joukowski airfoil. 
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Fic. 3. Contour for evaluation of leading-edge force. 
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Fic. 5. Elliptic cone in supersonic flow. 
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Fic. 6. 
line and over the airfoil chord between +1. An alter- 
nate treatment of the thin Joukowski airfoil is given in 
reference 1. 

As is evident from Fig. 2, the thin-airfoil theory shows 
a negative pressure over the forward three-fourths of 
the airfoil, changing with a constant gradient to a posi- 
tive pressure at the trailing edge. Such an increase of 
pressure toward the rear obviously leads to an upstream 
buoyancy, or negative drag, in contradiction to D’Alem- 
bert’s principle. A simple calculation shows that the 
buoyant force in this case (D = —volume X dp/dx) 


SCIENCES—MAY, 1950 


corresponds to a drag coefficient of 
Cp = —2r?? (5 


It will now be shown that this negative drag is exactly 
canceled by a positive contribution arising from the 
singularity in the flow field near the leading edge, pro. 
vided the integration of pressure over the airfoil syr. 
face is properly modified by a limiting process to jp. 
clude this singular term. To carry out this limiting 
process, we select, in preference to the airfoil surface 
itself, a streamline a short distance above the airfoj 
surface and, furthermore, allow the integration to ey. 
tend a short distance ahead of the airfoil nose. The 
drag is to be obtained from the product of the pressure 
and the slope of the streamline—that is, the product 
uw may be formed by noting that 


(u — tw)? = u? — w? — 2inw (6 


Considering a contour as illustrated in Fig. 3, the drag 
of a portion of the nose of the airfoil will be given by the 


integral 


D, = 1.P. of p sf (u — 1w)*dé = 


Tr 1 —¢ : 
I.P. of pV7# fi— = +2(Ve—1-—é&) 4+ ia 


(7 


The first term in the square brackets contains the sin- 
gularity. In the limit as the airfoil thickness ap- 
proaches zero and the contour approaches the real 
axis, this term contributes the value 


a ee Er 
Co, = 1? at - - d& = 2rt? (8) 
c (E — 1)(E4+ 1) 


over a small length of the path c just above the point 
—1. This value is equal and opposite the value ob- 
tained [Eq. (5)] from the approximate surface pressure 
distribution. The remaining terms in the integral (7 
need not be considered, since they correspond to this 
previously determined value. 

Further insight into the significance of this calcula- 
tion may be gained by an examination of the velocity 
field represented by the complex velocity function (1). 
As illustrated in Fig. 4, the velocity function shows that 
the horizontal velocity u approaches — © at the pout 
—1. Obviously, the leading edge of the airfoil will lie 
some distance to the left of this singular point, since the 
additive velocity « becomes equal and opposite to the 
stream velocity V before the point £ = —1 is reached. 
The point at which u = —V is the stagnation point 
and marks the position at which the stream divides to 
form the upper and lower surfaces of the airfoil. Be 
hind the stagnation point the velocities given by the 
thin-airfoil theory are, in fact, the velocities of am 
internal flow along the central source-sink distribution 
of the airfoil. It is evident that the velocity distribu- 
tion along the source-sink line does not represent a good 
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Pressure distribution over thin elliptic cone in supersonic 
flow. 


Fic. 7 
approximation to the actual surface velocity distribu- 
tion in the vicinity of the rounded nose. The exact 
shape and pressure distribution over the rounded nose 
are thus not determined by the thin-airfoil theory. 
The determination of the drag of the nose portion with- 
out knowledge of its éxact shape or pressure distribu- 
tion is possible, however, because of the invariant prop- 
erty of the contour integral. 

Similar considerations apply in the case of a lifting 
airfoil. 
proximation to a lifting airfoil of finite thickness, and 
the leading-edge force (a forward thrust in this ease) is 
determined from the singularity in the flow around the 
The plate in this case represents a 


Here, an inclined thin plate is used as an ap- 


edge of the plate. 
line of discontinuity, or a vortex sheet, along the 
interior of the airfoil. 

In either case the application of the theory to com- 
pressible flows is made on the assumption that the con- 
figuration of the flow exterior to the airfoil remains es- 
sentially similar to the corresponding incompressible 
potential flow. The contour of integration (Fig. 3) 
is further supposed to lie in a region where the additive 
velocities are small, so that the velocity u and the 
corresponding pressure can be corrected by the Prandtl- 
Glauert rule. 


ELLIPTIC CONE IN SUPERSONIC FLOW 


As an example of a supersonic flow, consider the case 
of the flattened elliptic cone which was previously con- 
sidered by Squire. Adopting a variation of Buse- 
mann’s method,* the velocity field may be described by 
an analytic function of the variable ¢ where 


a 


¢ = 2e/(1 + e?) (9) 
and 
e= (y + 22)/(xn + V x y* — 8°) (10) 


As pointed out in reference 3, the variable ¢ is the argu- 
ment of the general solution of Prandtl’s equation, of 
zero degree (cf. references 4 and 5 for the equivalent 


i 


solution of Laplace’s equation). The variable ¢ is an 
analytic transformation of ¢ that maps the interior of the 
Mach cone onto the whole complex plane, transforming 
the Mach cone itself into the positive and negative 
branches of the real axis outside the points +1. The 
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variable ¢ approaches the undistorted space coordinate 
(y + iz)/x near the plane of the airfoil inside the points 
+1 and thus facilitates the choice of functions to 
represent various airfoil boundary conditions. Figs. 
5 and 6 illustrate these relations. 

As in the case of the Joukowski airfoil, the elliptic 
cone requires a singularity of the order —1/2 in the 
expression for the vertical velocity w to represent the 
rounded leading edge. A trial shows that the simple 
function* 


w= —(Vt/2)(m Vm? — *) (11) 


satisfies the boundary condition for an elliptic cone of 
maximum thickness xt having its leading edges along 
the lines y = mx. 

The horizontal perturbation velocity u, which deter- 
mines the pressure distribution (Ap = —puV), is 
found from w with the aid of the condition for irrota- 

that is, Ou/O0z = Ow/Ox. In terms of 


the variable ¢ this relation is equivalent to 


tional motion 


du/dt = i(¢/V1 — ¢2) (dw dé) 


or 


. 


“= if (¢(/V 1 — £°)(dw/dt)dé (12) 


where the integration begins at a point of zero disturb- 
ance. 
Integration for the ordinate of the surface yields 


* w * w = 
z= / pax = -y f ve aq = 3 V x? — (y/m)? 


(13) 


The constant ¢ is seen to be the maximum thickness of 
the section at x = 1.0. 

In the present example, the sweep angle is assumed 
to be greater than the Mach angle so that the elliptic 
cone lies entirely within the Mach cone. The veloci- 
ties uw and w then vanish on the Mach cone, and u 
can be determined by integrating Eq. (12) from —1 to ¢ 


Vt m _ +2 
ihess oe Fy \ . s+ F(®,k) — E(®,k) 
21— m- c2 — m? 
(14) 
where 
? l1-—¢ / ~ 
® = are sin \ . b= W1 — ow? 
1 — m°* 


and F and E are the elliptic integrals. A plot of the real 
solution for u, which is proportional to the pressure 
distribution, is shown in Fig. 7. It is interesting to 
note that the elliptical cross section leads to a constant 


* The proper value of the function is selected with the aid of 
the condition that w is to be discontinuous across the real axis be- 
tween +m and positive on the upper side 








surface pressure distribution. This result has _pre- 
viously been obtained by Squire. 

As in the preceding example, the drag is determined 
by integrating the product —puw along a stream sur- 
face just above the airfoil surface. Eqs. (11) and (14) 
for u and w actually contain two (i.e., the real and 
imaginary) solutions. Since the real solutions are the 
ones used here, it is necessary to separate the product 
of the real solutions from the complex product ww. 
With the aid of Eq. (12), it can be shown that the re- 
quired product of the real solutions is equal to one-half 
the real part of the complex product in a limited region 
near the nose of the airfoil where the term ¢/V 1 — °° 
can be considered constant. 

With this latter point in mind, the expression for the 


leading-edge force of a section at x = 1 can be writ- 
ten 
~ Vtm (V1l-& 
D, = R.P. p s x 
é 21-— mm V7 — m? 
Vt m 
: dé 


9 ' sea 9 
“Vm — ¢ 


seve - me V1 — ¢? 
= R.P.i? : J a a a 
4 L—m’* Sc (€ — mye + m) 


, m* : 
= p \ (15) 
2 {Vl — m 
or, in coefficient form, 
Co, = m(t?/4)(m/V 1 — m?) (16) 


for both edges. 

The contour for Eq. (15) crosses only a limited region 
near the airfoil nose and involves in the limit only the 
singular term ¢V 1 — ¢2/(¢2 — m?). Over the remainder 
of the interval the contour is equivalent to an evalua- 
tion of the surface pressure times the surface slope. 
Over this interval the pressure is constant, proportional 
to 


F(r/2, k) — E(r/2, k) 


and the resulting drag is easily computed from the pro- 


jected area of the cone—that is, 


—— T m — ; i 

partial Cp = — # - (F — E) (17) 
2 1— mm 

In this case, both the drag arising from the approxi- 

mate surface pressure distribution and the leading- 

edge force act in the same direction so that the total 

drag amounts to the sum of Eqs. (16) and (17). 


c ( t y| m* m$ 
‘D = 2 F —_ E S 
a NS V1 — m? - mt } — 


The quantity ¢/2m is the “thickness ratio” of the cross 
sections of the elliptic cone. Eq. (18) does not, of course, 
include the wake drag that would arise if the cone were 
cut off so as to have a blunt base. 
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Pp 2 
v =(VcosA) 
Fr =r 
I-(McosA)® 
Fic. 8. General expression for force per unit length normal to 


oblique leading edge. 


GENERAL EXPRESSION FOR LEADING-EDGE ForcE 


In the foregoing example, the flow field is conical, 
and the drag coefficient of each section throughout the 
length of the conical body is the same. In more general 
examples such as that described by Lighthill,* however, 
this simplifying feature will not appear, and the co- 
efficient of the leading-edge force will vary from point 
to point along the length of the edge. Furthermore, it 
will not be possible, in general, to represent the entire 
geometry of the flow field by a function of a complex 
variable. However, if the body is sufficiently thin and 
flat and presents no abrupt changes in the radius of 
curvature of the leading edge, the flow field in the im- 
mediate neighborhood of a point on the leading edge 
may be considered to have a cylindrical, or two-dimen- 
sional, form. In such cases the localized region of the 
flow field can be represented by a function of a complex 
variable, and the contour integration for the leading- 
edge force can be performed. Fig. 8 illustrates the more 
general relation obtained when this calculation is ap- 
plied to an elliptical leading edge with nose radius 1. 
It is seen that the force normal to the edge at large 
angles of sweep is approximately equal to that de- 
veloped by the impact pressure of the normal stream 
component acting on the developed area of the circular 
cylinder having the radius of the leading edge. 


REFERENCES 


' Allen, H. Julian, General Theory of Airfoil Sections Having 
Arbitrary Shape or Pressure Distribution, N.A.C.A. Report No 
833, 1945. 

2 Lighthill, M. J., Methods for Predicting Phenomena in the 
High Speed Flow of Gases, Journal of the Aeronautical Sciences, 
Vol. 16, No. 2, pp. 69-83, February, 1949. 

3 Jones, Robert T., The Use of Conical and Cylindrical Field 
in Supersonic Wing Theory, from compilation of papers presented 
at the N.A.C.A.-University Conference on Aerodynamics, June 
21-23, 1948, Langley Field, Va., pp. 341-353. (Durand Re- 
printing Committee, California Institute of Technology, Pasa- 
dena 4, Calif.) 

‘Hobson, E. W., The Theory of Spherical and Ellipsoidal 
Harmonics, pp. 164, 165; The University Press, Cambridge, 
1931. 

5 Bateman, H., Partial Differential Equations of Mathematica! 
Physics, pp. 8356-358; Dover Publications, New York, 1944. 





of 

occ 
for 
all 

wel 
air 
des 
pla 
sch 


cut 

con 
giv 
gra 
gen 
blu 
api 
sta! 
ver 
e.g 

util 


illu 
tior 
be | 
vic 
noz 
phe 


adj 


vel 
inc 
det 
ten 
col 


tio 
tie: 


tee! 
enc 
bee 


Co 


nee 


al to 


the 
lex 
ng- 
ore 
up 


rge 


1m 
lar 


‘0. 


{ne 


Cs, 


al 
re, 








A Sharp-Focusing Schlieren System’ 


ARTHUR KANTROWITZ{ anp ROBERT L. TRIMPIT 
Cornell University 


SUMMARY 


The conventional schlieren system, employing a parallel beam 
of light, is unable to distinguish between the density gradients 
occurring at various positions along the light path. The image 
formed by such a system thus represents the integrated effects of 
all the density gradients across the region under investigation, as 
well as the other sections of the light path such as glass walls and 
air outside the region of investigation. In many cases it would be 
desirable to obtain an image in which density gradients in a given 
plane determine the image obtained—i.e., a sharp-focusing 
schlieren system. : 

A schlieren system with multiple sources and corresponding 
cutoffs has this sharp-focusing property. Each source-cutoff 
combination produces an independent schlieren image, and, for a 
given screen position, only the shadows produced by density 
gradients at a single plane in the flow superpose exactly. In 
general, for planes out of focus, offset superposition of the images 
blurs out the effects of density gradients not in the focal plane by 
a process similar to the focusing of ordinary lenses. In its present 
state of development, the sharp-focusing schlieren system is not 
very effective near opaque walls parallel to the optical axis 
e.g., for boundary-layer studies—since only grazing rays can be 
utilized in this case 

The performance of a model sharp-focusing schlieren system is 
illustrated by a series of photographs of a jet and a flame posi 
tioned along the optical axis. It is shown that photographs can 
be taken of the jet undisturbed by the presence of the flame and 
vice versa. A series of photographs of a normal shock wave in a 
nozzle are shown and compared with a conventional schlieren 
photograph. 

Some of the more important problems in construction and 


adjustment are discussed. 


CONVENTIONAL SCHLIEREN SYSTEM 


is A CONVENTIONAL SCHLIEREN SYSTEM the light devi- 
ations produced by optical density gradients are con- 
verted to shadows. The shadow production is nearly 
independent of the position along the optical axis of the 
density gradient (Fig. 1), and thus final local light in- 
tensity depends on the sum of the density gradients en- 
countered by a ray in traversing a flow pattern. 
Therefore, as is well known experimentally, a conven- 
system has no focusing proper- 


tional schlieren 


ties. 
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PRINCIPLE OF OPERATION OF THE FOCUSING SCHLIEREN 
SYSTEM 


The focusing properties of a common lens are depend- 
ent on the divergent rays of light from the object. The 
convergence of these rays by the lens then forms the 
image. At locations other than the conjugate focal 
plane** to the object, the various rays originating from 
a point on the object do not superimpose exactly and 
the image formed is blurred; while at the conjugate 
focal plane the convergent rays superpose exactly 
(neglecting aberrations) so that a sharp image is formed. 
A similar superposition of images formed by divergent 


** That plane in which the image for a given object is located 
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Fic. 1. Sketch illustrating that, for a conventional schlieren 
system, the light intensity changes at the screen are independent 
of the position of the density gradient producing the shadows 
(A focusing lens is usually used with conventional schlieren 
systems and is useful for bringing solid objects into focus but does 
not focus the schliere.) 
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Fic. 2. Focusing effect with multiple sources. For density 


gradients occurring at a focal plane conjugate to the plane of the 
screen (as illustrated) the shadows produced by the schlieren 
system superpose exactly. For planes other than the focal 
plane, the shadows produced superpose with offset, thus blurring 
the image. This action is similar to image formation by ordinary 


lenses. 
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Fic. 3. Optical components of the sharp-focusing schlieren 
system used in this work. The sharpness of the focusing prop- 
erties can be judged by the angle e between the outermost di- 
verging rays. 





Fics. 4a-4d. Photographs of a high-pressure nitrogen jet 
5/¢, in. thick by °/,,in. high. (a, top left) Jet located at the focal 
plane. (b, top right) Jet located !/2 in. away from tbe focal plane. 
(c, bottom left) Jet located 1 in. away from the focal plane. (d, 
bottom right) Jet located 1!/2 in. away from the focal plane. 


beams of “‘parallel light” is the principle of the focusing 
schlieren. 

The sharp-focusing schlieren (see Fig. 2) employs 
multiple sources and slits with each source-slit pair a 
conventional schlieren system in itself. Although the 
beam of light originating from a particular slit may be 
considered as ‘‘parallel light’’—1.e., the rays comprising 
the beam are parallel—the beam itself is oblique to the 
optic axis of the system as it traverses the tunnel. Con- 
sequently, the light intensity changes that are caused by 
density gradients in the focal plane conjugate to the 
screen will superimpose exactly and thereby form an 
image, but the light intensity changes due to density 


gradients in other planes will superpose offset from each 
other and will be blurred out. It is still true that each 
beam has encountered various disturbances in travers. 
ing the flow and that its total deflection, which deter. 
mines the final light intensity, depends upon an jp. 
tegrated value of all these gradients. However, only 
the disturbances at the conjugate focal plane are com- 
mon to all the rays that arrive at a given point on the 
screen, and, hence, these disturbances alone determine 
the resultant sharp image. Fig. 2 illustrates this 
principle for two sources, and Fig. 3 shows schematically 
the actual optical units employed in our f« ycusing 
schlieren system. Mention should be made of the fact 
that, since the actual deflection of light beams produced 
by air density gradients is extremely small (its influence 
on the light intensity is attained only by means of the 
sensitive schlieren effect), one need not be materially 
concerned with distortion of the images by density 
gradients not in the focal plane. 


PHOTOGRAPHS OBTAINED WITH THE FOCUSING 
SCHLIEREN SYSTEM 


Figs. + show the focusing effect as a °/¢4- by * jg-in. 
nitrogen jet is moved away from the focal plane of the 
system. This jet was then placed 1!/s in. from an 
acetylene burner, and a glass sheet was placed between 
the two as shown in Fig. 5. The photographs resulting 
from this experiment appear in Figs. 6. From Fig. (ia 
(flame off) and Fig. 6b (flame on), it is evident that the 
presence of the flame does not interfere seriously with 
the observation of the jet when the focal plane is located 
in the jet. In Fig. 6c the converse is apparent, since 
the jet is almost invisible when the focal plane is at the 
flame. 

Figs. 7 show schlieren photographs of the shock wave 
in the flow through a 2- by 2!/2-in. nozzle at a Mach 
Number of approximately 1.3. In Fig. 7a, taken with 





Fic. 5. Photograph of the experimental apparatus with a 
setup used for photographing a flame through a jet and vice 
versa (see Fig. 6). The camera is off the photograph to the 
right, and light source and condensing system are seex: at the left. 
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A SHARP-FOCUSING 


the focal plane at the center of the channel, the entropy 
discontinuity extending downstream from the shock 
intersections is visible, thus demonstrating the sensi- 
tivity of the system. Note the absence in Fig. 7a of the 
mottled appearance (indicating turbulence) of the flow 
behind the shock wave which is usually present in con- 
ventional schlieren photographs (Fig. 8). 

Fig. 7b shows the resultant image when the focal 
plane is located one-third the distance from the center to 
the wall. The ‘‘doubled’’ appearance of the shock 
wave is an indication that the shock is oblique to the 
optic axis of the schlieren system, in contrast to the 
center of the tunnel where it is parallel to the axis. The 
identification of such ‘“‘oblique’’ shocks was made by 
previously setting up the system with the optic axis not 
perpendicular to the tunnel walls and focusing on 
the shock at the center of the tunnel where, by 
symmetry, the shock Was normal to the walls. Thus, 
the shock curves into a “‘dished’’ shape across the 
channel. 

The shock formation with the focal plane !/y5 in. from 
the tunnel wall appears in Fig. 7c. There are no sharp 
density gradients noticeable in this photograph because 
of the interaction of the shock with the wall boundary 
layer. The mottled appearance of the picture is due to 
the wall boundary layer, and the development of tur- 
bulence after the shock is evident. The conventional 
schlieren photograph, Fig. 8, of the 
includes most of the aforementioned features, but 
they are superimposed so that it is impossible to dis- 
tinguish the various shock configurations across the 


same shock 


channel. 

Boundary-layer investigation with horizontal knife 
edge is not handled adequately with the focusing 
schlieren apparatus. This limitation is due to the 
opaque tunnel surface, which interferes with the passage 
of rays oblique to it, since beams inclined toward the 





Fics. 7a—7c. 


knife edge. 





The oblique shocks originate at a wall that is just below the photographs. 
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Fics. 6a (left) and 6b (right). Photographs, with the setup 
shown in Fig. 5, are taken with the focal plane at the jet and with 
the flame off (a) and flame on (b). 





Photograph, with setup shown in Fig. 5, is taken 


Fic. 6¢ 
with the focal plane at the flame with the jet on 


surface will be reflected and those directed away from it 
will pass through the boundary layer only near the side 
walls. A mild focusing effect for the boundary layer 
may be obtained by the use of a few (probably even one 
would be sufficient) long horizontal slits. The elimina- 
tion, by masking, of the source slits that yield the re- 
flected rays was used in the experiments performed to 


date. 





2 


Sharp-focusing photographs of a shock wave, Mach Number 1.3 in a nozzle (2 by 2'/, in.) with a horizontal 


(a, left) is focused at the center of 


the channel, and the sharp normal shock wave, as well as the entropy discontinuities following the shock intersections, can be seen. 


(b, center) is focused one-third of the distance from the center toward the wall. 
(c, right) is focused !/,. in. from the wall. 


shocks are oblique to the optic axis at this stage. 


The doubling of the shock wave indicates that the 
The mottled appearance usually associ- 


ated with nozzle flow photographs is seen to be due to flow at the wall boundary layer rather than in the center of the tunnel for the 


region covered by these photographs. 
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Fic. 8. Photograph of the same shock pattern taken with 
a conventional schlieren system. In this photograph most of the 
effects of the previous three photographs can be seen to be super- 
posed, thus making it difficult to distinguish the plane in the 
flow where the various density gradients occur. 


SUGGESTIONS FOR CONSTRUCTION AND ADJUSTMENT 
TECHNIQUES 


Since one must effectively align multiple schlieren 
systems simultaneously, a high degree of accuracy in 
both the manufacture and operation of the system is 
imperative. The focusing schlieren system employed 
was designed both mechanically and optically to 
minimize the number of required adjustments and 
still retain the desired freedom of motion for precise 
alignment. 


(1) Optical Elements 


After earlier versions of the system employing simple 
lenses for collimation were unsuccessful, it was deter- 
mined that anastigmatic collimating lenses (Fig. 3) 
were a necessity to minimize the distortion of the source 
slits furthest from the optic axis. In the present sys- 
tem, good quality wide-angle photographic lenses are 
satisfactorily used for this purpose. 


A lens not usually found in the conventional schlieren 
setup is required for the focusing apparatus. Since an 
enlarged image of the spark source is needed to illumi- 
nate the entire field of source slits, a condensing system 
is used to form this image. However, the outer diverg- 
ing rays resulting from this condensing system will not 
strike the first collimating lens unless deflected toward 
the optic axis. This deviation is caused by the addi- 
tional achromatic lens as shown in Fig. 3. 

The second achromatic lens, located after the cutoff 
plates, accomplishes the necessary focusing of the 
images on the screen. 


(2) Source and Cutoff Plates 


The heart of the sharp-focusing schlieren system jg 
the multiple source and cutoff plates. The distance of 
the outer active (illuminated) source slits from the optic 
axis determines the maximum obliquity of the rays to 
this axis. Since this obliquity is necessary for sharp 
focusing, it is desirable to have as large a distance as 
possible between the outermost slits. (The limit for 
this distance is fixed by the accuracy of manufacturing 
anastigmatic lenses to overcome the distortion of the 
outermost slits.) Also, as has been previously dis. 
cussed, the effective blurring out of nonfocused dis. 
turbances is dependent on the number of superimposed 
images, and, consequently, a large number of slits js 
requisite. <A sufficient number of slits can be attained 
without undue difficulty. 

The source plate employed consisted of 82 slits, of 
which about 50 were illuminated, and even a smaller 
number might have been sufficient providing the same 
maximum obliquity was retained. These slits were 
1'/: in. long, approximately 0.004 in. wide, and spaced 
about 0.036 in. apart. These are not critical dimensions 
so long as the cutoff matches accurately with the source 
plate (see below), and appreciable variations from these 
dimensions occurred with no noticeable disadvantage. 
The. method of manufacture to ensure accurate match- 
ing of the source and cutoff plates will be described 
after the requirements for the plates are presented. 

After several unsuccessful attempts with mechanical 
methods, it was decided that a photographic reduction 
procedure (onto glass plates) could most easily produce 
the source and cutoff plates with the desired accuracy, 
while satisfying the following criteria. The first re- 
striction imposed was that the emulsion of each plate 
face the test section, so that deviation of the light in 
passing through the glass of the plates could be ignored. 


(Continued on page 319) 
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MOTION FOR 
CUT-OFF 


KNIFE EDGES 

Fic. 9. Schematic illustration of source and cutoff plate con- 

figurations. To ensure accurate simultaneous cutting off of all 

of the multiple schlieren systems, a technique of producing the 

source and cutoff plates must be used which ensures that like 

cross-hatched edges are of identical contour. This is accomp- 
lished by photographic techniques discussed in the text. 
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RIEF REPORTS of investigations in the aeronautical sciences 

and discussions of papers published in the JOURNAL are 

presented in this special department. The publication will be com 

pleted approximately 8 weeks after receipt of the material. The 

Editorial Committee does not hold itself responsible for the opinions 

expressed by the correspondents. Contributions should not exceed 
800 words in length. 


On Wing-Body Interference at Supersonic 
Speeds 


George K. Morikawa 

Missile Aerodynamics Sectidn, Hughes Aircraft Company, Culver City, 
Calif. 

February 8, 1950 


N A RECENT PAPER, ! Ferrari has attempted to correct and mod- 
| ify his results arising from his study of wing-body interfer- 
ence in linearized supersonic flow.?_ In the more recent paper are 
included experimental results of pressure distribution measure- 
These data (in reduced form) are compared in Fig. 1* 

The body 
i.e., the in- 


ments. 
to a theoretical result based on a different analysis.* 
considered is an infinitely long circular cylinder 

fluence of the nose is neglected; the spanwise upwash distribu- 


* Ferrari’s notation is used 


tion yields slightly lower values than those produced by Ferrari's 
body in the vicinity of the wing. 

Some remarks concerning the present analysis, and, in par 
ticular, the theoretical spanwise pressure distribution shown in 
the figure may be appropriate. The analysis in the region outside 
the downstream Mach cone emanating from the leading-edge- 
body junction (called the ‘‘body’’ Mach coné) is exact within the 
assumptions of linearized flow theory; specifically, in the figure, 
the calculation is exact for y > 1.404. (The wing-tip correction 
is not shown, but this can be calculated by well-known methods 
for planar systems.) Within the ‘‘body’’ Mach cone, the analysis 
is markedly more difficult, but an approximate solution is obtain- 
able, as exhibited by the theoretical spanwise pressure distribu- 
tion in the figure for 1 < y < 1.404. It appears unfortunate that 
pressure measurements were not made at a larger chordwise dis- 
tance from the leading edge, since the region of particular inter- 
est—i.e., within the ‘‘body’’ Mach cone—is so small here. 

Several interesting observations may be made from Fig. 1: (1) 
The usual nonlinear behavior with angle of attack, 8, is shown by 
the spread of the four experimental curves; (2) the relative in- 
effectiveness of the expansion side of lifting surfaces is empha- 
sized; (3) for y > 1.404, the theoretical curve nearly coincides 
with the distribution given by the smaller angle of attack (8 = 
—4°) on the compression side; and (4) within the ‘tbody’? Mach 
cone, even the compression side of the wing appears to be less ef- 
fective than the theory predicts for small 8. It is hoped that 
more experimental work will be available in the near future to 
compare with the available theoretical work on nonplanar sys- 


tems. 
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1 Ferrari, C., Interference Between Wing and Body at Supersonic Speeds 
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Stabilization of Shock Waves in Channel Flow 


Arthur Kantrowitz and Stanford E. Neice* 

Graduate School of Aeronautical Engineering, Cornell University, 
Ithaca, N.Y 

February 22, 1950 


ie CONVERGING-DIVERGING SUPERSONIC DIFFUSERS, the normal 
shock has a stable equilibrium position in the diverging portion 
of the diffuser. 
always present in practice, move the shock away from its equilib- 
rium position, expansion pulses moving it downstream and com- 
pression pulses moving it upstream. When a compression pulse 
is large enough to force the normal shock past its unstable equi- 
librium position in the converging portion of the diffuser, then 
the shock will continue to move upstream, thus ‘‘unstarting’’ the 
diffuser. For maximum efficiency, it is, of course, desirable to 
reduce the shock intensity to a minimum 
consistent with the requirement 


i.e., move the shock 
as close as possible to the throat 
that the flow must not be unstarted by the compression pulses 
ordinarily present. To achieve this aim, the use of a long throat 
was proposed in reference 2. Various experimental confirma- 
tions of the effectiveness of long throats in improving the ef- 
It should 
be remarked that these efficiency increases are obtained in spite 
of the increase of skin friction caused by the long throat. 

It occurred to us that a similar stabilization of a shock wave 
could be accomplished by the use of a surge chamber connected 
to the throat of a supersonic diffuser through slits in the walls, 


ficiency of supersonic diffusers have been obtained. 


as is shown in the upper half of Fig. 1. 

When the diffuser is operating normally with the normal shock 
in its equilibrium position, the velocity at the throat will be 
supersonic, and the pressure in the surge chambers will correspond 
to the low static pressure of that region. If a compression pulse 
moves the shock beyond the throat, the reservoir openings will be 
exposed to the high-pressure air behind the normal shock. The 
resultant airflow into the surge chambers will act to remove part 
of the additional mass from the compression pulse. The result 
will be to decrease the upstream velocity of the shock. If the 
initial pulse is not too large, it will be brought to rest before it 
passes the forward unstable equilibrium position in the converging 
part of the channel, thus preventing the unstarting process. The 
function of the surge chamber in removing the excess mass is 
needed only occasionally when a strong compression pulse comes 
from the rear. After the passage of the pulse and the return of 
the shock to its normal position, the surge chambers are open to 
the low supersonic channel pressure and so return to their initial 
condition. The energy necessary to accomplish this return to the 
initial pressure comes from the channel flow and is taken from 
the flow only when the shock is downstream from the slits. Thus, 
it reduces the available energy of the supersonic stream and tem- 
porarily moves the shock slightly upstream, reducing the strength 
and losses of the shock. Thus, the work required to return the 
surge chambers to their initial pressure comes from energy ordi 
narily wasted in the supersonic diffuser. 

* Graduate Student in Aeronautical Engineering This note is based on 
a thesis submitted in partial fulfillment of the requirements for the degree 
M.Aero.E The authors are grateful to the Office of Naval Research for 


partial support of this research 
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Disturbances coming from the rear, which are 


SCIENCES—MAY, 1950 





THROAT 
VOLUME —s RESERVOIR 
VOLUME 

SK 
ae, ee a ae ; 














vo i t 
FY KN Wwwy 7 
\ = 
SHOCK STABLE PULSE 
POSITION — ~UNSTABLE TUBE 
VELOCITY 
> PROFILE 
O 
O 
a 
LJ 
> 
CHANNEL STATION 
Fic. 1. 
We define: 
V, = volume of the channel between the stable and unstable 
equilibrium position of the normal shock 
V. = volume of the surge chambers 
p1 = average density of the air in the throat corresponding 


to the shock in its stable equilibrium position, super- 
sonic flow 

p2 = average density of the air in the throat after the shock 
has passed the slits and approached the unstable 
equilibrium position 

pe, = final chamber density after shock has passed the 
slits and approached the unstable equilibrium posi- 


tion 

Pe, = initial chamber density, shock in stable equilibrium posi- 
tion 

T) = stagnation temperature; it is assumed that the temper- 


ature of the surge chamber is always stagnation tem- 
perature 
7, = stream temperature before shock 
The excess mass in the throat caused by the passage of the 
shock is approximately V;(p2 — p:), and an air mass of this order 
must be removed by the chamber to produce a significant im- 


provement. 
The increase in mass in the chamber is [since p.. ~ pi(71/T 
equal to V.(p., — pe,). 
Equating these quantities and dividing by pi, we have approx! 
mately 
(T1/T 0) Ve[(pe,/pe,) — 1] = Vil(pe/pi) — 1] 
If we designate p.. = p-(1 + m), where m is the fraction 


2 1 
increment of p.. (or the fractional increment of the pressurt 


then 


Ve = (10/71) (1/n) Vi [(p2/p1) — 1] 


This expression gives the volume required in terms of the frac 
tional increment of the initial density, assuming absorption ol 4 
mass equal to the excess mass in a pulse that would just unstart 
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It will be desirable to keep » small to 


, conventional diffuser. 
avoid large openings into the diffuser passage with attendant 


Josses. 

For practical ground installations, the reservoirs required to 
keep » small will not present any difficulty and can easily be 
made so large as to absorb the required mass without any appre- 
ciable pressure increase. For air-borne installations, on the 
other hand, the space available for surge chambers will be re- 
stricted. It may be desirable in such cases to use a larger slit 
area than in the ground installations, so that the required mass 
can be absorbed even though the reservoir pressure rises appre 
ciably. The cross-sectional area of the required slits will be com 
parable but smaller than the second throat area. 

In order to analyze the required slit dimensions, the method of 
characteristics as applied to nonstationary one-dimensional flow 
has been extended to include the effect of mass removal. With 
the application of this extended method of characteristics, it is 
possible to plot the position of the shock on a time-displaceme.at 
diagram and to observe the effects of mass removal. A series of 
such calculations, made assuming an infinite reservoir, indicated 
that openings 40 per cent of the throat area were required to 
double the ability of the channel to absorb compression pulses. 
These slits would permit the channel to absorb, without unstart- 
ing, compression pulses of twice the maximum pulse area that 


could be absorbed without the slits. 


EXPERIMENTAL RESULTS 


Experimental work on this device was done in a small blow- 
down tunnel with a supersonic diffuser throat 1.21 by lin. The 
maximum Mach Number before diffusion was 1.65. The inves- 
tigations were undertaken to see whether an increase in the maxi- 
mum pulse that could be absorbed by a supersonic diffuser could 
be increased by the use of surge chambers. To this end, slits 
40 per cent of the throat area, as indicated by the theoretical 
analysis of reference 1, were placed about one throat diameter 
downstream from the minimum section. This location permitted 
the surge chambers to absorb mass over a longer interval of 
time, since the shock moved past the slits earlier, thus increasing 
the slit effectiveness. 

Observations were made, first, under disturbances normally 
present in the channel. Visual observations and schlieren pic- 
tures with continuous illumination, made with the normal shock 
just downstream from the slit, showed that the shock motion was 
greatly reduced by the presence of the surge chamber, as was 
expected 

A pulse tube, much like a miniature shock tube, with an ex 
tremely short compression chamber, was constructed to produce 
a triangular pulse and was placed near the point where the dif 
Thus, a part of the pulse 


fuser discharged into the atmosphere. 
The strength 


from the pulse tube propagated up the diffuser. 
of the pulse generated could be varied by varying the pressure in 
the compression chamber. Starting with the channel shock at a 
definite known position and with the slits closed, it was found 
that the diffuser could be unstarted if a gage pressure of 25 lbs 
per sq.in. or greater was used in the compression chamber. When 
the slits were opened, it required a pressure of 65 Ibs. per sq.in. 
in the compression chamber to unstart the diffuser. Calculations 
show that the pulse area generated by the shock tube in the sec- 
ond case was roughly twice that generated in the first case. This 
is in agreement (indeed in fortuitously close agreement) with our 
theoretical expectations. 

It may be concluded that considerable stabilization of the 
position of a normal shock wave can be obtained by the use of 
surge chambers. Comparative measurements of maximum ef- 
ficiency should be made in the pressure of disturbances normally 
met in practice to evaluate the relative merits of this and other 


schemes for stabilizing supersonic diffusers. 
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An Observation Concerning the Rankine- 
Hugoniot Relations 


Jack Lorell 

Research Enginzer, Jet Propulsion Laboratory, California Institute of 
Technology, Pasadena, Calif. 

February 27, 1950 


kes RANKINE-HUGONIOT RELATIONS imply that the entropy 
increment, AS, across a normal shock wave is positive, going 
downstream. The proof of this fact seems to be glossed over in 
many of the reference sources on supersonic flow and shock 
theory. For example, Liepmann and Puckett, on page 41 of their 
text,! include a proof valid for only weak shocks. 

Perhaps the following proof of this entropy relation is of in- 
terest. Let subscript (); refer to upstream, let subscript ()2 
refer to downstream conditions, and let the symbols y, M, S, 
and C, carry Then, the following 


equation applies across a normal shock wave: 


their usual meaning. 


y-1 2 -—1\'% 
. (: + ue) FA = 
AS 2 : oe. +7 1 
= In “4 =In@ 
’ a SS 
(1) 


[ef. Eq. (4.7) in Liepmann and Puckett]. The quantity ¢ in 
Eq. (1) equals 1 for 14, = 1. Furthermore, if M,? — 1 = m > 
0 (and the difference 7,2 — 1 need not be small), then d @/dm ex- 
ists and is positive. It follows that, for all values of Mj, > 1, 
AS > 0. 

The assertion concerning d¢/dm becomes obvious on evaluating 


this quantity: 


2(7 — 1)m? 


: ional. 2y 
(y +1 {1+ m (1 + m} (1 + m) 
y¥+ 1 y¥+1 


1 do = 


gdm 


Thus, d¢/dm > 0 provided only m > Oand y > 1. 
It may be of interest, also, to point out that the minimum 
Mach Number after a normal shock can be easily obtained 


i.e., from the relation 


1 y-1 1 1\ 
M? — M,? — — =(-+. (3) 
2y 24 2° 24 


[ef. Eq. (4.4) in Liepmann and Puckett]. Thus, 


M,? > (y — 1)/2y7 
since the right side of Eq. (3) is always positive. 
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The Influence of Some Parameters in the 
Flutter Problem 


G. lsakson and T. F. O'Brien 
Aero-Elastic and Structures Research Laboratory, Massachusetts Insti- 
tute of Technology, Cambridge, Mass 


Received February 2, 1950 


B™ AND ARNOLD! HAVE SHOWN that flutter of an elastically 
supported airfoil is theoretically possible for the limiting 
case of zero air speed when the nodal point of a natural mode 
of vibration of the airfoil is located three-quarters of the 
chord aft of the leading edge. They also show that the exist- 
ence of zero air-speed flutter requires that structural and air 
friction be absent and, for this condition, that the flutter speed 
increases as the nodal point is located at increasing distances from 
the three-quarter chord point. 
when structural friction is allowed, the flutter speed reaches a 


However, they show that, even 


minimum when the node is close to the three-quarter chord point. 

These results suggest that a study of the relationship between 
node-point location and the inertia and elastic characteristics of 
the airfoil, and particularly of the characteristics that determine 
a node at the three-quarter chord point, might provide informa- 
tion that would be useful in designing so as to avoid flutter. This 
is the purpose of the present brief investigation, which, it is 
hoped, will stimulate further investigation along similar lines. 

A simple model, as shown in Fig. 1, is first assumed. It con- 
sists of a rigid airfoil connected by a tension compression spring 
and a torsion spring to a rigid fuselage which is free to become 
displaced vertically but is not free to pitch. 
nection of the two springs to the airfoil represents its elastic cen- 
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The frequencies, w-, and @,, of the normal modes of free vibra- 
tion of the system, or the so-called coupled modes, may be deter- 
mined by solution of the equation, 

‘ 2 2 2 2 
ee CID *S25)*"1+ 
Wh p? + a? Wp Wh my \p? + a? 
(3) (+) 
1+ =0 (1 
Wh my 

where w; and «» are the frequencies of the restrained torsion and 
bending modes, respectively, when the fuselage is held stationary; 
my» and my are the masses of the airfoil and fuselage, respectively; 
p is the pitching radius of gyration of the airfoil about its center 
of gravity as a fraction of the chord; and a is the distance of the 
airfoil center of gravity aft of the elastic center as a fraction of 
the chord 

When the natural frequencies have been determined, the cor 
responding mode shapes may be expressed in the form, 


0c p ‘ . . 
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where @ is the angle of pitch, positive nose up and ; is the vertical 


displacement of the airfoil, positive down. C is the chord length. 
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It is easily seen that the distance of the node point aft of the 
elastic center, /, as a fraction of the chord, is given by the equa- 
tion, 


l = —(6c/z)— 3 


Eqs., (1), (2), and (3) thus show that the node location depends 
upon the parameters p, a, mw/my, w:/w», and the chordwise loca- 
tion of the elastic center. When the node is located at the three 


quarter chord we have 
1 = (3/4) -—e (4 


where e¢ is the location of the elastic center aft of the leading edge, 
as a fraction of the chord. Introducing Eq. (4) into Eq. (3), we 
obtain the following condition for the node to be at the three- 
quarter chord: 


e = (3/4) + (0c/s)—! 5 

When my/my, p, and w;/wy) are held constant, Eq. (5) expresses 
a relationship between the elastic center location and the center- 
of-gravity location. This relationship is shown plotted in Fig 
2 for a number of values of w;/w, and for representative values of 
Mw/my and p. The value chosen for mw/my; was 0.18, which is 
roughly equivalent to a value of 1 for the ratio of wing to fuselage 
weight for an actual airplane with tapered wing. This equiva- 
lence is based on a representative section at 70 per cent of the 
wing semispan and dynamic similarity of the model and airplane 
for vibration in the fundamental bending mode. The value 
chosen for p was 0.25. 

When a = 0 and when the frequency ratio w/w, = [1 + 
(mw/my) |'/2, the frequencies of the two coupled modes are equal 
For this special case the two coupled modes are not uniquely de- 
termined and may consist of any combination of bending and 
torsion. Consequently, the elastic axis locations giving a node 
at the three-quarter chord point are indeterminate. It may be 
noted that for w/w, < [1 + (mw/my) |'/2, the coupled mode that 
is predominantly torsion has a lower natural frequency than the 
coupled mode that is predominantly bending. 

Fig. 2 provides some interesting indications. It shows, for 
example, that, if low-speed flutter is to be avoided, the ratio 
«w/e», should be as large as possible—that is, the torsional stiffness 
of the wing should be great in comparison with the bending stiff- 
ness. 


the elastic center, but not too far aft, and that the elastic center 


It shows also that the center of gravity should be aft of 
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READERS’ 


itself should be as far forward as possible. These conclusions 
are not new. They have been stated before as the result of more 
elaborate analyses and an extensive accumulation of data. It is 
interesting, nevertheless, that they can be reached so simply from 
the criterion of Biot and Arnold. 
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A Sharp-Focusing Schlieren System 


(Continued from page 314) 


Then, it was also necessary, as shown in Fig. 9, to satisfy 
the requirement that each cutoff edge approach the 
matching edge of the source image—i.e., the contour of 
the cutoff edge be identical to that of the opposite edge 
of the source image for which it was acting as a knife 
edge. 

The source and cutoff plates were made by photo- 
graphing, against a white background, a piece of trans- 
parent plastic on which were fixed various strips of 
black scotch tape. Scotch tape '/16 in. wide, spaced at 
\/,.in. intervals, was photographed, and this negative 
was used as the source plate. 

The plastic sheet was then turned over, additional 
strips of scotch tape '/, in. wide were affixed so that 
one edge of the wider tape coincided roughly with the 
centerline of the '/,-in. tape, and another photograph 
was taken of the acetate in its new position. The nega- 
tive obtained was the positive of the required cutoff 
plate and, consequently, had to be ‘“‘reversed’”’ by proc- 
essing or by exposing another negative held behind it in 
parallel light. 

Another possible method of manufacturing the plates 
again obtains the source plate by photographing '/j¢-in. 
scotch tape. Then, without turning the plastic sheet 
over, the '/,-in. tape is attached and a second photo- 
graph taken. This negative is used to expose the final 
cutoff negative by a contact exposure process in which 
the emulsions of the two negatives must be touching. 
If other processes of producing the source and cutoff 
plates are used, care must be taken to ensure that the 
corresponding edges (those edges cross-hatched in Fig. 
9) match, since a severe loss in performance will result 
when nonidentical edges approach each other during 
“cutoff.’’ 

Kodalith Orthochromatic and Kodak Process Pan- 
chromatic plates have both been used successfully for 
the sources and cutoffs, with the former yielding a 
slightly higher quality plate. Precise focusing of the 
camera was accomplished by viewing, via a microscope 
and a hole cut in the rear of the plate holder, the image 
formed on the emulsion of an exposed plate. 


(3) Mechanical Design 

The mechanical design of the system was determined 
to allow: (a) the alignment and adjustment of the 
optical system in itself; and (b) the movement of the 
entire optical system so that complete coverage of the 
tunnel was possible even with the restricted field (1 in. 
minimum diameter) of the lenses. Vernier adjustment 
was provided for the movement of the source and cutoff 
plates parallel to the optic axis. Rotation about an 
axis parallel to the optic axis and motion perpendicular 
to the optic axis of the cutoff plate were also governed 
by verniers. The complete optical system was mounted 
on rails and adjustable pedestals to permit coverage of 
the entire flow field. 


(4) Adjustment Techniques 


The highlight of the final adjustment technique, 
assuming that all optical units have been previously 
aligned, will now be briefly described. The source and 
cutoff plates are placed at approximately the focal 
distance from their respective collimating lenses. 
Then the longitudinal vernier attachments are manip- 
ulated until the image of the source slits exactly super- 
pose on the corresponding cutoffs. A rough longitudinal 
adjustment may be obtained by eliminating the dark 
bands that appear on the cutoff plates when the source 
image is of a different size from the cutoff negative. 
During the above manipulation, rotational adjustment 
must also be made to correct nonparallel source-cutoff 
conditions. It has been found that the system is 
correctly adjusted when the image formed on a screen 
darkens uniformly, while the cutoff plate appears to 
darken radially (i.e., darkening progresses evenly to- 
ward the center from all directions) as the knife edges 
are moved to “‘‘cutoff.’’ The focal plane of the system 
is determined by placing an object, such as a fuzzy 
piece of string, at the desired plane and then moving the 
image screen until this object is in focus. The best 
photographs have been obtained with the light only 
slightly ‘‘cut-off,’’ since under this condition the effects 
due to irregularities in the plates are minimized. 
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